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ABSTRACT 


Requirements for submarine periscope depth operations have been increased by 
integration with carrier battle groups, littoral operations, and contributions to joint surveillance. 
Improved periscope depth performance is therefore imperative. Submarine control personnel 
rely on a large number of analog gauges and indications. An integrated digital display system 
could enhance the ergonomics of the human control interface and display additional 
parameters. This thesis investigates the required feedbacks for robust automatic depth control 
at periscope depth, and thus indirectly determines the additional parameters desired for an 


integrated display. 


A model of vertical plane submarine dynamics is coupled with first and second order 
wave force solutions for a particular submarine hull form. Sliding mode control and several 
schemes of state feedback are used for automatic control. Head and beam seas at sea states 
three and four are investigated. The automatic control effectiveness provides insight into the 
indications used by the ship's control party in operations at periscope depth. One possible 
display system is proposed, with several additional enhancements to improve ship's safety, 
reduce operator fatigue, and enable accurate reconstruction of the events leading to a loss of 


depth control. 
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I. INTRODUCTION 


E GENERAL 

The need for attack submarines to operate at periscope depth has been increased by 
integration with carrier battle groups, operations in the shallow littoral, and contributions to 
joint surveillance. 

Operating at periscope depth beneath a seaway, a submarine is in an unstable 
condition. Аз Ше free surface 1s approached, the seaway forces increase, trying to pull the 
submarine to the surface. To counter these forces, the ship’s ballast is changed and control 
surfaces are used. Because of the seaway's stochastic nature, manual operation for long periods 
at periscope depth taxes the ship's control party. 

Operators must remain aware of the environmental conditions. If the sea becomes 
quiescent, the submarine will sink out. If the sea suction forces are greater than the ballast and 
planes authority, the submarine will broach the free surface increasing detection risk by several 
orders of magnitude. Other events, like temperature or salinity changes, can also have major 
effects on reliable depth keeping. Contributing to the environmental issues, the need to use 
minimum speed for a given sea state to control the detectable mast feather reduces the 
available planes authority, and increases the difficulty of depth control. 

However, the current submarine force is not optimized for these operations. One 
inexpensive area for improvement is the display system for the ship's control party. Modern 
digital display systems offer ergonomic improvements over current gauges and readouts. 

Given a requirement to conduct submarine ship control manually, a fundamental 
question is that of how to display the state of the ship to the operators. Aside from the 
obvious indications like ship’s pitch angle, depth, and control surface positions what else would 
be useful? Candidates include the net force acting on the ship, accelerations, and various time 
averaged values. Implied in this is that a nontraditional means of display will be used to show 
these parameters, so that the operators will not have to rely on a number of gauges or meters, 
with averaging of results only available only by the calibrated eye. 

An intelligent assistant to the ship’s control party would show items of current 
concern, and issue alerts based on an operator programmable doctrine. Issues like mast 
exposure, ship’s relationship to the bottom, and trim state could be shown in an intuitive, 


logical manner. 


Current evolutions and other items relating to the tacucal employment could be included as 


required. 


B. AIM OF THIS STUDY 

Although the ship’s control party currently relies on a small number of indications, the 
ability to sense “by the seat of the pants” cannot be discounted. This thesis investigates 
required feedbacks for robust automatic depth control at periscope depth, and thus indirectly 
evaluates the additional indications to be added to an integrated display. 

This approach assumes that the best ship's control partes already use system states for 


control which are not explicitly displayed. 


C. THESIS OUTLINE 

Chapter II contains the development of the deeply submerged submarine dynamics 
model. Chapter III gives the development and source of the wave forces used to simulate 
operations at periscope depth. In Chapter IV, optimization studies are performed for nine 
different cases of state feedback control. This gives a feeling for the quality of depth control 
achievable by the use of different levels of sensors. Chapter V explores the use of sliding 
mode control for periscope depth operations. In Chapter VI, current ships control technology 
Is reviewed and an integrated display is proposed. Conclusions and recommendations are 


given in Chapter VII. 


II. SUBMARINE DYNAMICS MODEL 


A. INTRODUCTION 

When a submarine is deeply submerged, many of its maneuvering characteristics can 
be determined from application of Morison’s equation to model test data. A series of trials, 
often done with a planar motion mechanism (PMM), give the damping and inertia coefficients 
for small maneuvers in each of the six degrees of freedom. This method is not without limits. 
For trials done in the horizontal and vertical planes only, nonlinear cross coupling effects are 
ignored. The hydrodynamic coefficients work poorly for prediction of high speed maneuvers 
and control surface casualties. Here the large crossflow velocities, vortex hull interaction, and 
flow separation all have effects which are not predicted by the hydrodynamic coefficients. It is 
possible, however, to include some of these effects as additional nonlinear terms. 

As the submarine approaches the free surface, several complexities are introduced into 
the hydrodynamic coefficient approach. First, the inertia terms change as an acceleration will 
no longer act upon an effectively infinite region. Second, an inviscid form of damping exists 
near the free surface. This comes about from the generation of waves by the body, and 
depends on the body depth and character of motion. Finally, the interaction between the 
incident waves and the submarine introduces added forces and moments. These effects 
combine to make designing for periscope depth vexing for engineers and operating at 
periscope depth an art for the ship’s crew. 

The approach in this thesis will be to first establish a dynamics model appropriate for a 
deeply submerged submarine at low to moderate speeds. The forces and moments resulting 
from the seaway will then be superimposed on this model to provide a reasonable 


approximation to the submarine motion beneath waves. 


B. DEEPLY SUBMERGED EQUATIONS OF MOTION 

К Definition of coordinate system and states 

The coordinate system defined in Figure 1 will be used. The origin of the global 
coordinate system is fixed at the ocean surface. The z axis is positive downward, towards the 
ocean bottom. The x axis is positive in the direction of intended submarine motion. The body 


fixed coordinates are rotated from the global coordinates by the angle 0. Body fixed velocities 


w (heave) , u (surge), and q (pitch) are shown. The control surface deflections, ó, (bow planes) 


and 6, (stern planes) are also defined. 


Ocean Surface 





Figure 1. Coordinate System Definition 


2 Hydrodynamic coefficients review 

For a deeply submerged submarine, small motions can be analyzed using the concept 
of hydrodynamic coefficients. These represent a Taylor series expansion of the functional 
relationship between body movements and the resulting fluid forces. For example, given the 
deeply submerged body in Figure 2 undergoing pure heave, resulting body forces can be 


expressed in the following manner: 


и 


Figure 2. Submerged body in pure heave 


M = M,w+ M a ww] M, w ob 


Z = Z,w-t Z4 ww] Zi vw (2) 


This method is extended to the six degrees of freedom of the body, and done for 
velocity and acceleration components of the movement. This includes representations of 
added mass, viscous drag, and square law drag. 

В: Vertical plane equations of motion 

By using this system of notation, and applying Newton's second law to the body fixed 
coordinates, and transforming to global coordinates, the equations of pitch and heave may be 
obtained in the vertical plane. The general case is quite complex, having centers of mass and 
buoyancy that are separate from each other and the coordinate system origin. This, along with 
cross coupled hydrodynamic coefficients, results in a nonlinear, coupled set of differential 
equations. 

These equations of pitch and heave may be simplified considerably by several 
reasonable assumptions. Assuming that the submarine motion is constrained to the vertical 


plane, the equations of motion for heave and pitch are (Smith, Crane, and Summey (1978)): 


ту — uq - xcd – 269 ]= (5) 
Zo qct Wa ug i 2 w 


ES u*(Zg 6, +Z50,) 


I,q -m[|xç (w—uq)— zç (u +wq)] = (4) 
M ,q * M, wet M,uq * М шо 
*u^ (M5 Ó, * М5б,) 
— (xo mg — xg B)cos(0) 
— (zo mg — zgB)sin(0) 


It is apparent that Equations (3) and (4) are nonlinear, coupled differential equations in 
» and 4 and 4. To reduce this coupling, terms involving the derivatives of w and 4 can be 


collected, resulting in a mass matrix. 
М = 


The mass matrix can be readily inverted: 


^ m (6) 


w 


"VES 


By applying Equation (6), the cross coupling of terms in w and q can be removed 
from Equations (3) and (4). To allow the introduction of external forces and moments, the 
system was augmented by force and moment disturbances acting at the origin of the body 
fixed coordinates. They were multiplied by the cosine of the pitch angle for conversion to the 
body fixed coordinate system. These disturbances can be used to input external effects, such 
as changes in trim and wave forces. By further assuming that the center of buoyancy 1s at the 
body fixed coordinate system origin, the center of mass is directly below, and that the forward 


speed z is constant, the equations of motion can be reduced to the following: 


W = a, uw +a,,ug +a), Sin(0) * bj u^ 6, - bu? 6, * F, cos(0) +e,,q" +€,,qw (7) 
q = a,yuw + d55uq + a», sin(O) +b; u ó, + Бри, + M,cos(0) + eq T €34qW (8) 
Ó = q (9) 

Z 2 wcos(0) — usin(0) (10) 

x = wsin(@) + ucos(@) (11) 


where: 


2 a LIU М+М, 
" (m-Z,X1, - Mj)-Z,M, 


ааа cae Metu ш 
8 (т-220,-М,)-2,М; 


E MAN Zoe te ned 
5 MZA p M,)-Z,M 


E] | 
e MO Tm )М. 
% (т-2.Х1,-М,)-2,М; 


ауу = MES 
EE CZ SME 


He Q ¿r.o CNE 
"^ (m-Z,)u, - M)-Z,M, 


_ Uy Mis + 25M, 
“(те 2 Пре Мун М; 


p a MsZ + Ма(т- 2.) 
"^ (m-Z,X1,- M)-Z,M, 


Т” ақыт аа а 
7 (т-2,Х1,- М,)-2,М, 


b M,Z; % М;(т-2,) 

2 (т-2,Х1,-М,)-2,М, 
a (1, - M,)zgm 

" (m-Z,1,- M))-Z,M, 


Z¿2¿M 


OS SAPE 
(m- Z, X1, - M,)- Z,M, 


C» М „гст 
"^ (m-Z,),- M) - Z,M, 


= uz zm 
“ (т-2,Х1,-М))-2,М, 


Equations (7) through (11) are the governing equations of motion for this thesis. It is 
of note that the disturbance force and moment terms represent accelerations due to the 
disturbances. To provide ease of use, the equations of motion were implemented in the 
SIMULINK? model shown in Figure 3. This building block approach was very effective for 


conducting studies on the effectiveness of different types of controllers. 


Block info and requirements 


. Net local 
disturbance vector 


Mux [integrator Output 
State Vector 


Control Surfaces Control 
Matrix 


nonlinear 





Figure 3. SIMULINK® model of vertical plane submarine dynamics 


For control design, it is convenient to use a linear state space representation of the 


system. This allows the use of a variety of controller design tools including pole placement and 


linear quadratic regulator algorithms. Equations (7) through (11) can be linearized about a level 


flight condition. This results in the linear state space representation: 


W = a, uw ta, ug +4,,0 +b, u"d, +b,,u76, + F, (12) 
d - ajuw * a uq * a440  b4u4 6, but 6, М, (13) 
0-q (14) 

z=w-ud (15) 

x=w@+u (16) 


Equations (12) through (15) can be rewritten in matrix form. This form of the linear 
submarine vertical plane dynamics equations will be used for controller design. For controller 
design, Equation (16) was excluded from the matrix form. Because of the constant forward 


speed y assumption, there was no direct means of control for x. 


w| [аш apu аз О|у| |Бји“ Бри“ Е, (17) 
4| |ази ауи а, 0|4(| |5,5 byu? f | |Ma 
өзі ут” НИК НЕ 0 
2 l 0 и 0/2 0 0 0 


OF EXTENSION TO VERTICAL PLANE PATHKEEPING 
Equations (7) through (10) and the corresponding SIMULINK® model are linearized 
around a constant commanded depth, or level flight. They can be extended to a two 


dimensional pathkeeping simulation by a coordinate transformation. After coordinate rotation 


by an angle D (positive in the same direction as 8), the resulting system is: 


w 2 a, uw aug * at, sin(0')+bu26 + F, (19) 
à 2 a, uw * aug + ai, sin(0’)+b,u25+ M, (19) 
Ө%-4 (20) 
¿"= wcos(0”) — u sin(0”) (21) 
e CE омой (22) 
where: 
9'-0- (23) 
x'= -zsin(B) + xcos() (24) 
z' 2 zcos() * xsin(B) (25) 


| (26) 


a,, =a,,cos(B) 


E (27) 


a,, =a,,cos(f) 


Е, = F, +a,,cos(0”)sin(P) (28) 


M, = M, +a,, cos(0”)sin( ) (29) 


If the expected angular deviation from the planned path is small, Equations (28) and 
(29) can be simplified by assuming that cos(0^) is equal to one. Then the rotated equation set, 
Equations (18) through (22), 1s identical in form to Equations (7) through (11). 

Equations (23) through (29) allow any vertical plane path consisting of a series of 


straight line segments to be simulated one segment at a ume. 


D. THE DARPA SUBOFF 

1. Background 

For the purpose of this work, it was desired to have a vertical plane model of 
submarine dynamics which would give a similar response to a modern fast attack nuclear 
submarine (SSN). Several sets of unclassified hydrodynamic coefficients were available, these 
being for the swimmer delivery vehicle (SDV) detailed in Smith, Crane, and Summey (1978) 
and for the DARPA SUBOFF model detailed in Roddy (1990). 

The SDV had a very complete set of hydrodynamic coefficients which have been used 
in a large number of Autonomous Underwater Vehicle (AUV) research projects. Among these 
is the Naval Postgraduate School (NPS) AUV sliding mode controller, Hawkinson (1990). 
Despite these advantages, the SDV hydrodynamic coefficients were not used because the wing 
like hull of the SDV bore little resemblance to an axisymetric submarine hull. 

The SUBOFF hydrodynamic coefficients detailed in Roddy (1990) lacked some of the 
cross coupling coefficients. The documentation also lacked details on the models metacentric 
height. Because the SUBOFF represented a submarine hull form and most of the vertical 
plane coefficients and parameters were available, it was chosen as the model for this thesis. 

2. SUBOFF known parameters and coefficients 

The SUBOFF was developed to allow comparison between flow field predictions and 
model test data (Roddy, 1990). The available coefficients were based on planar motion 


mechanism tests conducted on the model. 
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Because the aim of the study was to examine full scale submarine motions, the model 
and its hydrodynamic coefficients were scaled to a length of 300 feet. After scaling, several 
parameters had to be modified or assumed to give control and response comparable to a 
modern fast attack submarine. The force coefficients of the stern planes was doubled to 
provide a more realistic level force. Bow planes were assumed to have one half the force and 
one quarter the moment authority of the stern planes. Finally, a metacentric height of one foot 
was assumed, as it provided a realistic point of stern planes reversal. The resulting parameters 


are shown in Table 1. 
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Table 1. SUBOFF Assumed and modified parameters 




















E. CONCLUDING REMARKS 

A simplified model of submarine vertical plane dynamics was derived. The coefficients 
for use in this model were obtained from the DARPA SUBOFF model, which is a 
representative axisymetric submarine hull form. The simplified nonlinear equations of motion 
were incorporated in a SIMULINK® model to allow easy integration with wave force models 


and different controllers. 
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Figure 4. DARPA SUBOFF model, Roddy (1990) 
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Ш. WAVE FORCE MODELING 


А. INTRODUCTION 

As a submarine operates near the free surface, it encounters complex forces which may 
cause unsatisfactory or unstable depth contro]. The lift and moment from incident waves 
increase in an exponential manner as the surface is approached. To maintain a desired depth, 
the ship’s ballast is adjusted to counteract steady forces. Control surfaces are used to counter 
dynamic changes. A small depth excursion or change 11 forces can overwhelm the planes and 
cause a loss of depth control. The consequences range from losing radio reception to 
compromising the ship’s mission. 

The effects of incident waves on a submerged body can be divided up in several 
categories. The largest, the first order forces act at the incident wave frequency. These forces 
move the submarine, but usually result in oscillations about a mean state. Second order forces, 
which are the result of wave diffraction and wave interaction, have several different frequency 
components. 

Wave diffraction of a single frequency wave results in a steady force and a varying 
force at twice the wave frequency. The double frequency force is typically neglected, as the 
large inertia of the submarine effectively filters it. Interactions of waves at different 
frequencies also results in forces. These consist of a component acting at the sum of the wave 
frequencies and a component acting at the difference of the wave frequencies. The sum 
frequency force 1s typically neglected, as it 1s also filtered by submarine’s inertia. The difference 
frequency component results in a slowly varying force on the submarine. 

The slowly varying forces are the principle cause of difficult periscope depth control 
(Ni, Zhang, and Dai, 1994). They are compensated for using control surfaces and occasional 
adjustments of trim. 

During the design phase, engineering decisions are made which will determine the 
ship's ability to remain at periscope depth. Of these, the most critical are the height of the sail 
and control surface sizes. Every foot added to the sail gives a deeper periscope depth. Larger 
planes improve the operator's ability to compensate for changes in suction forces. However, 
these improvements are not without cost. The sail and other appendages are a large fraction of 
the total drag, and can restrict the ship's top speed. Larger movable control surfaces can 


adversely affect the high speed casualty recoverability ( Jackson, 1992, p. 15-9). 


The goal of this thesis 1s not to provide new tools for the designer, but rather new 
means to enhance control for the operators of current submarines. Due to this focus, 


simplified means of modeling the wave forces for a few specific cases will be used. 


B. REVIEW OF LINEAR DEEP WATER WAVES 

The pertinent features of linear deep water waves will be reviewed to provide 
background for the following secuons. The coordinate system used for the examples is shown 
in Figure 5. For the examples in this section, it will be assumed that the submarine 1s oriented 
with the bow pointing into the page. Consistent with the global coordinate system from 
Chapter II, the distance from the surface to the submarine centerline is z. The submarine 


diameter 1s D. 


Wave speed,c 


Surface elevation tD 


shown at t=O 





Figure 5. Coordinate Definition for plane progressive wave, adapted from Sarpkaya and 
Isaacson (1981, p. 151) 


For a wave of wavelength L, a wave number, &, can be defined. 


ae (30) 
E 


Assuming that fluid is incompressible and inviscid Laplace’s equation can be applied. 


It is thus desired to find a solution to: 


ә ә 2 
o^ a 


To this, the boundary conditions at the free surface, and the bottom must be applied: 


16 


= =0 ,at z=d (no flow through ocean bottom) eo 
2 + II + ° =0 , at z=7 (zero velocity normal to ocean surface) (9) 
X 
ab 1|(986у (офу (ол) 
Bet Dike Күт еле ЫШ 
For small amplitude waves in deep water, the following solution can be obtained 
(adapted from Sarpkaya and Isaacson, 1981, p. 159): 

ф = — e^ sintor) (35) 

2 о? =g (36) 

k k 
: D 
P (37) 
2л 

П = —cos(@t) (38) 

= -2 ek со$(ах) (52) 

C = oet sin( or) (40) 

> H 1 

£ = о? E соѕ( 001) ш 

Е= -Zet біл(ал) _ 

| Н 45 

tc SPICE cos(@r) ee) 

2 (44) 


H ={7 _. 
~ e = Sinton ) 
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where 7] 1s the distance from the surface to the average level (z 2 0), ois the angular 
frequency of the incident wave, ¢ 1s the displacement of a particle in the x direction, and ё 


is the displacement of a particle in the z direction. 


A key parameter in oscillating flows is the Keulegan-Carpenter number: 


K = О. (45) 
D 


where U, is the average velocity across the characteristic dimension D. 


жан 


By taking the average of the velocity given in Equation (40) , and substitution into 


Equation (45), the expression for the Keulegan-Carpenter can be reduced to the following: 


(46) 


Equation (46) 15 the Keulegan-Carpenter number based on the cross flow velocity of the 


undisturbed wave at the same depth as the centerline of the submarine hull. 


C. WAVE FORCE REGIMES 

There are different regimes of interaction between a submerged body and a wave field. 
Broadly, they can be broken into several areas. Inertial interaction, where the body acts like a 
particle in the wave field. Wave diffraction, where the bodies influence upon the wave field 1s 
accounted for. Finally, there are flow separation (viscous) effects. The relative importance of 
each of these effects can be determined by examining the relationship the body size to the 


wave parameters. (Sarpkaya and Isaacson, 1981, pp. 381-386) 
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Figure 6. Wave force regimes (Sarpkaya and Isaacson, 1981, pg. 385) 


To estimate the significant effects for a typical SSN, a typical operaung condition is 
assumed. For a 300 foot submarine with a 35 foot diameter, a typical periscope Operating 
depth would be about 50 feet from the centerline of the ship to the free surface. Using average 
values for sea states three and four and assuming deep water compared to the wavelength, the 


following quantities were calculated at a depth of 50 feet: 






Table 2. Estimated Wave Loading Parameters 






The Diameter/Wavelength (D/L) ratios and the Keulegan-Carpenter numbers of 
Table 2 show that for the sea states of interest, wave diffraction much more significant than 
viscous forces. It can be concluded that an inviscid analysis should give good results for the 
wave forces. However, this 1s only rigorous for an unappended hull, as the control surfaces 


and sail on an actual submarine will experience viscous effects. 
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D. SOLUTION FROM SLENDER BODY THEORY 

Wave force solutions for several specific cases were generated for the SUBOFF by the 
SSBN Security Department of the Johns Hopkins University Applied Physics Laboratory. A 
slender body solution with some three dimensional correcuons was used. The specific method 
used for the generation of the ftrst order motions and second order forces is detailed by O'Dea 
and Barr (1976, pp. 7-25). 

A seaway approximation consisting of a small number of regular waves was used to 
model sea states three and four. For each sea state, the resulting data were separated into two 
categories. The effects of the first order forces were given in terms of body motions. The 
effects of the steady second order forces and the difference interaction forces were provided in 
pounds force. 

1. Seaway model 

A random seaway can be represented by the superposition of a large number of regular 
waves. The seaway was approximated by superimposing z regular waves. The frequency and 
height of these waves was determined using the Bretschneider spectrum. It gives the spectral 
density in terms of the significant wave height, H, , and the peak frequency, o,. 

5H? d A | > 


S(Q —— — 
AT 


To model sea state three, a significant wave height of three feet was used, with a central 


frequency of 0.836 radians per second. This results in the following spectrum: 
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Figure 7. Example Sea State three spectrum 


Figure 7 gives a statistical picture of the seaway, but is not immediately useful for ume 
domain simulation. One way to obtain a time history is to represent this stationary process as 


a the sum of a series of sine waves: 


48 
n(t)= >, A, $11( 00.1 + а, ) (48) 


i-l 
Where A, is the amplitude of the 7^ wave , and ais its randomly chosen phase angle. 


If the number of sine waves 15 reasonable large, and the frequencies and amplitudes of 
each component are chosen to achieve the same energy as the section of spectrum it 
represents, Equation (48) will give a good representation of the ocean surface. 

The method chosen was to divide the spectra into n segments of equal areas. This 


results in n sine waves all with equal amplitudes. Integration of Equation (48) yields: 


ЕИ 


a e 
) с 


Because the spectrum extends to infinity, 1t was chosen truncate the spectrum at a 
point where the area was a fraction C of the total area. The amount of area to be represented 
by each sine wave 15 equal to its mean square value. So the amplitude of each sine wave is 


equal to the square root of the area it represents times the square root of two. 


Рини та (50) 


- 5 


ELE 


Equation (48) can be integrated up to some frequency w,, which represents the 


frequency at which the spectral area is equal to iC/n times the total area. 
(51) 


[5 4 ETT 
А (оао 16 n 


Solving Equation (51) for o, yields: 
(52) 


Because the spectral level is insignificant below «w equal to 0.6o, , the frequency if the 


first segment was determined as follows: 


_ (060, +01) (53) 


w, 2 


The remainder of the frequencies were determined by taking the midpoint of the 
frequencies at either side of the area segment. 
Q, +0 (54) 


0 ==, for i=210 n 
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Figure 8 illustrates the method used, approximating the spectrum with sinusoids. 


Nineteen equal area sections are divided, with the center frequency of each segment marked 
with a circle. 
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Figure 8. Spectra area division and mean frequencies 
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Figure 9 shows the ocean surface which results from the use of this method for the 
case of sea state three, peak frequency of 0.862 radians per second. Nineteen sinusoids were 


used to approximate the spectra, and the phase angles were randomly chosen. 
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Figure 9. Sea surface approximation for sea state three using nineteen sinusoids 


2. First order forces 
The first order wave effects were provided in the form of submarine motions. They 
were given as a series of phasors, the real part of the summation representing the actual 


perturbation caused by the first order wave forces. 


z(t) = AA (22 


Q(t) = S Decir ren (56) 
Because the first order motions were provided for a specific depth, it was required to 
correct Equations (55) and (56) for depth. The first order motions roughly correspond to the 


particle motions given by Equation (42), so an exponential decay was used to derive the 


following correction factor: 
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p^ (5) 





Applicaton of Equation (57) to Equations (55) and (56) results 1n: 
2(t) = У’ ge Marsa )-k G-z) (58) 


j=} ! 


0(t) = po = os C (59) 


The displacements given by Equations (58) and (59) are not suitable for inclusion in 
the submarine equations of motion. For this, an acceleration is required. Differentiating twice 


with respect to time results in: 


45 n 2 I: =k (>= 60 
ša) = - 2... O Z e i(í ta, )—-k (2-2,) ( ) 


é(t) = ES д.2. е-аи+а,)-к (гг) (61) 
iz} I f 


Equations (60) and (61) were incorporated as force and moment disturbances in the 
equations of motion found in Chapter II. To test the validity of this approach, an open loop 
simulation was performed using the accelerations from Equations (60) and (61) for one sea 
state and heading. Figure 10 shows the results of this simulation, as well as the expected first 
order motions. The upper curve shows the expected first order motions, and the lower curve 
shows the results of integrating Equations (7) through (11) with the accelerations from 
Equations (60) and (61). Although there was some drifting motion, the character of motion 
and the approximate amplitude of each cycle of motion very close. The drifüng motion was a 


result lack of the lack of open loop depth stability, which is characteristic of submarines. 
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Figure 10. Submarine response to first order accelerations, and expected response 


3. Second order forces 
For a particular depth and wave time history, the second order forces were given in the 


following form: 


A p ence rm) (62) 
[=] у=! 


м) = pus 2 мо „а“ еђна +a,)} (63) 
i=} g=) 7 


Z(t) represents the force acting at the body fixed coordinate system in the z direction 
and M(t)is the moment acting about the y axis. It should be noted that Equations (63) and 


(65) include the slowly varying forces (i  j) and the steady forces (i = j). 
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It can be determined from analysis of the changes of second order forces with respect 
to depth given by Crook (1994, pp. 61,62) that the steady forces with the following exponential 


decay factor: 


e ^* (64) 


The slowly varying order wave forces vary with depth according to the sum of the wave 


numbers: 


ook thy 2 (65) 


SOS 


Application of Equations (64) and (65) to Equations (62) and (63) results in: 


(66) 


дозе, Xr rene 
= gs 


М (1) = E pee ye -аа јичана,у- н, уга.) (67) 


The real portion of Equations (66) and (67) represents the steady and slowly varying 


second order wave forces acting on the submarine, with correction for depth. 


4. Inclusion of wave forces in equations of motion 
The first order accelerations and second order forces had to be combined to form the 
force and moment disturbance accelerations for use in the deeply submerged equations of 


motion (Equations (7) and (8)). 


ЕА PA P (68) 
= =| +M 
M „ (1) и O(t) M (t) 
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Е. CONCLUDING REMARKS 

An elementary review of linear wave theory was presented. The case of interest, a 
submarine at periscope depth, was examined to determine the salient elements of its interaction 
with the incident waves. The parameters suggested that the major features of the incident 
wave effects on the submarine could be determined by using a potential analysis with inertial 
and diffraction forces accounted for. 

The Bretschneider spectrum was used to determine the spectral density functions of 
the sea states of interest. For the purpose of time domain simulation, the spectrum was 
approximated by the superposition of a number of regular waves with randomly chosen phase 
angles. 

The first order force transfer function and second order forces response amplitude 
operators were provided for the SUBOFF for a nominal speed and depth. Approximate depth 


scaling was introduced to allow use at depths other than nominal. 
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ГУ. STATE FEEDBACK CONTROL AT PERISCOPE DEPTH 


А. INTRODUCTION 

ile State feedback control 

One popular means of control is to feed back the system states after the application of 
linear gains. System response of linear systems subjected to this type of control is predictable, 
and a variety of tools are available for control law gain selection. 

The ship’s control party on a submarine with conventional indications does not have 
the full state of the ship to operate from. Although the actual instrumentation may vary 
somewhat, in general a few analog indications are used in conjunction with a digital depth 
indication. For this reason, various levels of partial state feedback were used to evaluate the 
effects of missing indications. 

The use of different state feedback schemes was felt to be appropriate to model human 
operators. The treatment of airplane pilots as a control law “has come to be recognized as a 
quasilinear element for random-appearing tracking tasks related to piloting. At the same time, 
the pilot retains spectacular nonlinear gain changing, mode switching, and goal seeking 
precognitive control capabilities as yet only partially explored.” (Graham and McRuer, 1991, p. 
1093) In this context, it was assumed submarine “pilots” could be treated in a similar fashion, 
with feedback from each operating state determined with linear gains. 

The use of a first order lag was considered to model the combined human and control 
surface response time. It was found that reasonable lag values (on the order of a half second) 
had minimal effect on the control response and corresponding submarine motions. Because of 
the computational expense, the control response time was neglected. 


State feedback control of the linear system 


Х = Ах + Ви (69) 


where: 


ел а а 
pe Control mathe 
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хє $ State Vector 
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u € &" control vector 


ean De expressed as: 


и = Кү (70) 


where: 


Keni (71) 


The system given by Equation (69) subject to the control law given by Equation (70) has the 
following closed loop dynamics matrix: 


А. = (А + ВК) (72) 


The eigenvalues of the closed loop dynamics matrix will be related to the system 
stability and responsiveness. In general, the real portion of the eigenvalues must be negative 


for system stability. Also the more negative the eigenvalues, the faster the system response. 


2: SUBOFF simulation parameters 

Wave force data was available for the SUBOFF for four different cases. These were 
sea states three and four with head and beam directions. All were valid at a speed of six knots 
and depths greater than fifty feet. 

At six knots, the linear state representation used for eigenvalue determination and 


control law design 1s: 


№] [-0.0179 37101 00196 0O[w] [-00628 -04009 F, (t) (73) 
q 0.0006 -0.0680 -0.0034 Olg 0.009 —.0027 д, М а (1) 
207 е о о |5 0 
Z ] 0 — 10.1269 Ol z 0 0 0 
where: 
Im (1) = Ж г wee 
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trim 


All simulations were performed using the nonlinear equations of submarine dive plane 


motion: 
№ = ануу t aj5uq * aj4sin(0) +b, ,u76, +bu ô, + F, cos(@) + e, ar + €,,qw (74) 
q =a, uw + a,,ug + a>; sin(@) + Би", + bu, + M ,cos(@) + eq +e qw (75) 
д7 (76 
z = wcos(@) — usin(@) (77) 
x 2 wsin(0) - ucos(0) (78) 


The simulations were performed using a commanded depth of 55 feet and using a zero 
error initial state vector. Commanded pitch angle, heave and pitch rate were all zero. The 
depth was chosen to provide a good representation of actual submarine periscope operating 
depth. 

9. State feedback implementation with SIMULINK® 

The state feedback controller was implemented in the SIMULINK® model shown in 
Figure 11. This block was designed to use an optional feedforward signal, and also to facilitate 
the use of integral control (Both feedforward and integral control are discussed later in this 
chapter). Deflection limits were placed on the control surfaces. Control surface rate limits 
were not included, but could be easily added. These limits are of interest because of the 
relationships between control surface rates, hydraulic plant size requirements and noise from 


control surface operations. 
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Figure 11. State feedback control block diagram 


A SIMULINK® model was developed to incorporate the submarine dynamics of 
Chapter II, the wave forces of Chapter III, and the state feedback control law. Also included 
was a logical means of adjusting the submarine’s trim. This was done by adding ballast in units 
of thousands of pounds at the center of buoyancy, and shifting ballast from the forward trim 
tank to the after trim tank in units of thousands of pounds. The details of the trim model are 


shown in Figure 12, while the overall model 1s shown in Figure 13. 
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Figure 12. SIMULINKS trim model 
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Е SIMULINK® state feedback control submarine model 


4. Integral control on depth 


To apply integral control, an additional state is introduced. Equations (74) through 
(77) are augmented by: 


2 = 27 Ссотталаеа (7 9) 


which is used to provide state feedback. This forces the steady state value of z to zero. In 
general, this approach is satisfactory as long as the control effort does not become saturated 


and the eigenvalues of the integral state are slower than the state which is being zeroed. 


5. Feedforward of wave forces 

Given the wave forces values, control effort can be directly applied to eliminate the 
average depth error. With a constant disturbance, a steady state value of the depth error can be 
determined (Appendix B). Using the linear equations of moton, the steady state depth error 


can be written as a linear combination of the net force and moment disturbances: 


“ss ~ <commanded = С, Fy + С, Ма (80) 
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To eliminate the depth error, it is required to apply the control effort it applied: 


K. A (81) 


К а 


Equations (80) and (81) can be combined to give a matrix gain relationship between the net 


disturbance and the feedforward: 


K $ Kia С С К _ C, Kj, С.К. p (82) 
м са м,| |ск, СОКА |М 
24 d |2525 2% 24 d 


The state feedback control law with feedforward is: 


Ó, 
09 ee 


sp 


(83) 


It has been suggested ( Musker, Loader, and Butcher, 1988), ( Ni, Zhang, and Dai, 
1994) that effective periscope depth control can be achieved by feeding forward the average 
second order wave forces. Because wave forces are a dynamic disturbance and the feedforward 
was calculated for a steady disturbance, a filter was employed to cut out the high frequency 
components of the wave forces. The filter employed was a first order Butterworth filter with a 


cut off frequency &,,. The cut off frequency was initially chosen as one radian/second. This 


was well below the maximum frequency wave force components ( around 2.2 radians/second). 
Figure 14 and Figure 15 show the effects of the first order Butterworth filter on the wave 
forces at a depth of 55 feet in sea state three. It is apparent that with a cutoff frequency of ten 
radians per second, the filtered forces and moments are very close to the unfiltered. At the 
lower cutoff frequency, 0.1 radians per second, the filtered forces and moments are much 
closer to the average values. 

To implement the feedforward control law, it was assumed that the net external force 
and moment were known quantities. Equation (82) was implemented in the SSMULINK® 


model shown in Figure 16 while the complete system model is shown in Figure 17. 
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Figure 14. Filtered wave forces for sea state three (head seas) 
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Figure 15. Filtered wave moments for sea state three (head seas) 
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Figure 17. SIMULINK? model of system with feedforward term 


6. Optimization algorithm and parameters 
One difficulty of using partial state feedback is that conventional pole placement or 
linear quadratic regulator algorithms can not be used to determine the gains. The gains in 


question were selected randomly, and gain combinations which gave stable eigenvalues were 
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simulated. Because of the clamping on the maximum planes angle, some gains which yielded 
stable eigenvalues resulted in unstable ship control. 

Randomly selected gains certainly provide less than optimum depthkeeping. Because 
of this, each feedback case was optimized to provide the best case for a particular sea state and 
commanded depth combination. In conjunction with the feedback optimization, the trim was 
optimized. 

The MATLAB® function CONSTR was used to perform the optimizations. 

CONSTR uses the Broyden-Fletcher-Goldfarb-Shanno variable metric method, and supports 
constrained optimization problems. To prevent the optimizer from selecting unstable systems, 
a constraint was placed on the eigenvalues. The real part of the eigenvalues was required to be 
less than a maximum value, usually -10~. 

The objective of the optimizations was to reduce the root mean square (RMS) value of 
the depth error. For the basic state feedback control, the average depth was expected to differ 
somewhat from the commanded depth of 55 feet. Because of this, the objective for these 
optimizations was to minimize the RMS value of the difference between the depth and the 
mean depth. 

Because the optimizations were performed without regard for minimizing control 
effort and or rates, large gains with attendant control chatter was expected. Although control 
chatter is not consistent with normal submarine operations, it was neglected to provide a clear 


basis of comparison between the differing levels of feedback. 


B. FEEDBACK OF DEPTH AND PITCH ANGLE 

T. Basic control 

An elementary level of ship control can be conducted with the stern and bow 
planesman, each operating to control one particular state. The logical approach to this is for 
the stern planesman to control the ship’s angle, and the bow planesman to control depth. This 
results in the following control law: 


| (84) 


communded 


М = И 
o и 0 O 0 Kia Q — Q commanded 
Z 


sp commanded 


ids e commanded 


D 


ló| « ó,... (85) 


After a stable set of random gains was determined, the controller was optimized to 
minimize the deviation from the average depth. The formal optimization statement 


(Vanderplaats, 1984, p. 9) is: 








Minimize: 
г (86) 
Ге 7 Zmean y dt 
UA e 
: * 
where: 
z = depth , determined by nonlinear simulation 
El 
| wat 
0 
¿mean ES , 
n 
H = Ballast added to center of buoyancy, thousands of pounds 
F = Ballast shifted from forward to aft, thousands of pounds 
Subject to: 


real(eigenvalues( A.)) S Eq, (87) 


Deviation from the mean value of depth, vice the commanded was used because of the 
expected average depth error. 

This approach was used for each of the four sea state cases. For sea state three (head 
seas), the optimized response is shown in Figure 18. The results of the four optimizations are 
shown in Table 3. For the RMS error and maximum error, the optimized values are given, 
along with their percentage of the initial values. 

In all cases, use of the optimization resulted in reduction of the mean square depth 


error (measured from the average depth). Reduction of the maximum error was also achieved. 
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-0.0074 + 0.20961 
-0.0074 - 0.20961 
-0.0356 + 0.11441 
-0.0356 - 0.11441 


-0.0074 + 0.20961 
-0.0074 - 0.20961 
-0.0356 + 0.11441 
-0.0356 - 0.11441 


-0.0074 + 0.20961 
-0.0074 - 0.20961 
-0.0356 + 0.11441 
-0.0356 - 0.11441 















Maximum Error (feet) 
Eigenvalues -0.0074 + 0.20961 
-0.0356 + 0.11441 

-0.0356 - 0.11441 


-0.0074 - 0.20961 
= Values pai la 

























RMS Error (feet) 0.4550 (49%) 0.7549 (82%) 0.657 (53%) 1.23 (95%) 
Maximum Error (feet) 1.533 (62%) 2.03 (82%) 2.54 (66%) 4.15 (87%) 


-0.0388 + 0.23921 -0.0419 + 0.14641 -0.0419 + 0.21781 -0.0010 + 0.21941 


-0.0388 - 0.23921 -0.0419 - 0.14641 -0.0419 - 0.21781 -0.0010 - 0.21941 
Table 3. Optimized pitch and depth control law results and performance 






























-0.0042 + 0.39111 -0.0010 + 0.23461 -0.0010 + 0.33971 -0.0419 + 0.13581 
-0.0042 - 0.39111 -0.0010 - 0.23461 -0.0010 - 0.33971 -0.0419 - 0.13581 
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Figure 18. Simulation with depth and pitch angle control in sea state three (head sea 
direction) 


2. Disturbance feedforward 
The pitch angle and depth feedback control can be implemented with a disturbance 


feedforward to correct average depth error. This results in the following control law: 


W — W commanded (8 8) 
дь, E 0 0 0 Ki Q — Q commanded H С, Ка С; Ki Е) 
О 10 0 K5 0 0 - mandes 0 0 Ма 
< — C commanded 
8] < б. (89) 


where Kiis given by Equation (82) and the force and moment disturbances are filtered. 
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After a stable set of random gains was determined, the controller was optimized to 


minimize the deviation from the average depth. The formal optimization statement 15: 





Т (90) 
| La C commanded y dt 
ык жс 
i 
Subject to: 
real(eigenvalues(A.)) < E nax (91) 


This approach was used for each of the four sea state cases. For sea state three (head 
seas), the optimized response is shown in Figure 19. The results of the four optimizations are 
shown in Table 4. For the RMS error and maximum error, the optimized values are given, 


along with their percentage of the initial values. 
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Initial Values 
Кэ 


0) (rad/sec) 


Eigenvalues -0.0074 + 0.20961 
-0.0074 - 0.20961 
-0.0356 + 0.11441 
-0.0356 - 0.11441 


МО ИУ +. LL ANI A S E s 7 
ERA EA ОРНА 
E TA S 2% CEST 
4 2 p е 
E AR 


Ku 
K> 


Eigenvalues -0.0337 + 0.33451 
-0.0337 - 0.33451 
-0.0092 + 0.60881 
-0.0092 - 0.60881 


-0.0074 + 0.20961 
-0.0074 - 0.20961 
-0.0356 + 0.11441 
-0.0356 - 0.11441 


-0.0354 + 0.60551 
-0.0354 - 0.60551 
-0.0076 + 1.04901 
-0.0076 - 1.04901 


-0.0074 + 0.20961 
-0.0074 - 0.20961 
-0.0356 + 0.11441 
-0.0356 - 0.11441 


-0.0322 + 0.86041 
-0.0322 - 0.86041 
-0.0107 + 1.62961 
-0.0107 - 1.62961 


-0.0074 + 0.20961 
-0.0074 - 0.20961 
-0.0356 + 0.11441 


-0.0250 + 0.56751 
-0.0250 - 0.56751 
-0.0179 - 1.601% 
-0.0179 - 1.60191 





Table 4. Optimized pitch and depth control law with disturbance feedforward results and 


performance 
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Figure 19. Simulation with depth and pitch angle control] with disturbance feedforward, sea 
state three (head seas) 


E. Integral control 

The feedback of depth and pitch angle can be augmented with integral control on 
depth to remove the average depth error. Since the bow planes are principally used for the 
control of depth, the integral control was applied to the bow planes only. This results in the 


following control law: 
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(92) 


W — Weommanded 


Q — Q commanded 
K 
15 0-6 
соттапаеа 


Sh 100 0 К, 
ба] [00 Ка 0 
% 2 m И титл 

21 


l| « ó... (95) 


After a stable set of random gains was determined, the controller was optimized to 


minimize the deviation from the commanded depth. The formal optimization statement is: 


Minimize: 





(94) 


£r 
24 
(2 Е £ commanded ) t 

0 


F(Kj4, Kis, K&,, H,F)— f 
f 


Subject to: 
(95) 


real(eigenvalues( ÀA.)) < E max 
This approach was used for each of the four sea state cases. For sea state three (head 


seas), the optimized response is shown in Figure 20. The results of the four optimizations are 


shown in Table 5. 
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^ 
щл 
к= 
— 
~J 
Un 
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щл 
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RTF OT pom 
SS Er 
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-0.0077 + 0.20881 -0.0077 + 0.20881 -0.0077 + 0.20881 | -0.0077 + 0.20881 


Figenvalues 
-0.0077 - 0.20881 -0.0077 - 0.20881 -0.0077 - 0.20881 -0.0077 - 0.20881 
-0.0318 + 0.11481 -0.0318 + 0.11481 | -0.0318 + 0.11481 | -0.0318 + 0.11481 
-0.0318 - 0.11481 -0.0318 - 0.11481 -0.0318 - 0.11481 -0.0318 - 0.11481 
-0.0070 -0.0070 

























Optimized Values 





MF LETT 


% 


— = оаа ааа 

== 
= 
Е 


Eigenvalues -0.0149 + 0.88241 -0.0134 + 0.52211 -0.0001 + 0.26151 | -0.0163 + 0.42961 
-0.0149 - 0.88241 -0.0134 - 0.52211 -0.0001- 0.26151 -0.0163 - 0.42961 
-0.0274 + 0.38871 -0.0286 + 0.24751 | -0.0421 + 0.16751 | -0.0258 + 0.17071 





Ki 
Kis 
Ka 


(RE FS PER 
4 д5» 
о РЕЯ 
. 
. 
. 


CE oue dA. 
Зе X k: LOFT, АХ dec d ERR. Т 
. . 
. 


S и 
5 жу; y А zE а 
» я 
. 
. 
. 


-0.0274 - 0.38871 -0.0286 - 0.24751 -0.0421 - 016751 -0.0258 - 0.17071 
-0.0012 -0.0020 -0.0015 -0.0015 





Table 5. Optimized pitch and depth integral control law results and performance 


45 


theta degrees 


500 1000 1500 2000 0 500 1000 1500 2000 


а mn nc 
; Цу 
ИИ 
x || mw 
ТЕ 11 ШИШ | 
|! | 2 INIM pJi | ны; 
; atun 
lee ia n spl rona fy 
АЈА A ТШ Ири 


ы ! || 


e as 
M d ua i: тї 


delta rad 


"iii 
m иу у 


ТЕ a, 
l 


EUIS luci vl 


| 
1 
| 
ШЕ у иши | 
| 


feet/sec degrees/sec 


500 1000 1500 2000 d 500 1000 1500 2000 
x feet x feet 





Figure 20. Simulation with depth, pitch angle, and integral control, sea state three (head 
seas) 


с. FULL STATE FEEDBACK WITH PARTIAL DISTRIBUTION 

1. Basic control 

The poor depth control of the previous section can be improved be adding to the 
number of states fed back. In keeping with previous logic, the bow planes will be controlled 
by the depth and heave, while the stern planes will be controlled by pitch and pitch rate. This 


results in the following control law: 
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na W commanded (96) 


Opp = Kj 0 0 Ki 4 — d commanded 
по к к. оо 


commanded 


>= 
< £ commanded 


5 = Sna (97) 


After a stable set of random gains was determined, the controller was optimized to 


minimize the deviation from the average depth. The formal optimization statement is: 


Minimize: 
(98) 


ЕКІ, Ку, Ку, Коз, Н, Е) = 





Subject to: 
real(eigenvalues(A.)) < E max (ОО) 


This approach was used for each of the four sea state cases. For sea state three (head 


seas), the optimized response 1s shown in Figure 21. The results of the four optimizations are 


shown in Table 6. 
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S5 


Sea 3/head 3/beam 4/head 4/beam 
State / Direction 


Initial Values 2 fe 
X% ; а 


i zd n NC CE M a A : 3%? И Б 2223 
KT 
22.5268 
24.9900 
0 
Mean Depth (feet) 55.06 55.33 5521 5579 
КІ 







НЛ (їй нав 


RMS Error (feet) 0.216 0.425 0.412 


Eigenvalues -0.6743 -0.6743 -0.6743 -0.6743 
-0.0413 + 0.35871 -0.0413 + 0.35871 | -0.0413 + 0.35871 | -0.0413 + 0.35871 
-0.0413 - 0.35871 -0.0413 - 0.35871 -0.0413 - 0.35871 -0.0413 - 0.35871 

-0.1789 ; -0.1789 





Optimized Values Е. 
SE ARAS А 
0 
0 


|77035 13.265 | > 2” 
1366.3 0 | 91224 91.49 o | 9644 
11633 o | 2087 51.15 о | 82.43 
89.26 0 3.813 0 21048| 0 15.33 0 


1.17 (147%) 0.99 (72%) 1.13 (90%) 3.46 (102%) 


Eigenvalues -3.3197 -0.1266 + 0.47091 -0.7286 -0.2213 + 0.88611 
-0.0695 + 1.25221 -0.1266 - 0.47091 -0.1458 + 0.47001 -0.2213 - 0.88611 

-0.0693 - 1.25221 -0.2619 + 0.01791 -0.1458 - 0.47001 -0.3797 + 0.48191 

-0.1879 -0.2619 - 0.01791 -0.1514 -0.3797 - 0.48191 





Table 6. Full state feedback (partial distribution) control law optimization results and 
performance 
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Figure 21. Simulation with full state partial distribution feedback control, sea state three 
(head seas) 


2, Disturbance feedforward 
As before, the basic control law can be modified to include a feedforward term to 


correct the average depth error. 


W — W commanded (1 00) 
МЕК 0 0 21 Q — Qcommanded 4^ Kia С, y [| 
de 0 K, Къз 0 US S 0 0 Ма 
€ — Z commanded 
[6] « ó,... (101) 
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After a stable set of random gains was determined, the controller was optimized to 


minimize the deviation from the average depth. The formal optimization statement 15: 





Minimize: 
7 (102) 
| (2 = геоттапаеа ) 4! 
PESE. Ke CH S РЁ) = | г, 
Subject to: 
real(eigenvalues(A,)) S Emax (103) 


This approach was used for each of the four sea state cases. For sea state three (head 
seas), the optimized response is shown in Figure 22. The results of the four optimizations are 


shown in Table 7. 
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12.5543 


Ye 


a o o y Се og 902 ү” 
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нз ASAS Po AA 


A AAA 


Exgenvalues 












22.5268 
24.9900 











-0.0413 + 0.35871 -0.0413 + 0.35871 | -0.0413 + 0.35871 | -0.0413 + 0.35871 
-0.0413 - 0.35871 -0.0413 - 0.35871 -0.0413 - 0.35871 -0.0413 - 0.35871 
-0.1789 


Optimized Values 


KT 25.81 0 15.73 0 12.74 0 16.87 
0 175.77 0 230.65 | 0 18.09 0 9 
0 26.1 0 53598 | 0 20.1 0 18 
3.5847 0 9.446 0 2.70 0 3.714 
BEEN. 


H/F (10? pounds) 7:17 799 198/8255 20.0/0.0 19.7/0.0 

Mean Depth (feet) 55.01 35 55.18 39:09 

RMS E:ror (feet) 0.0785 (21%) 0.3624 (68%) 0.5171 (94%) 1.25 (89%) 
Maximum Error (feet) 0.246 (19%) 1.07 (58%) 2.03 (118%) 3.26 (93%) 


Eigenvalues -0.1692 + 1.36011 -0.6730 -0.9159 
-0.1692 - 1.36011 -0.0429 + 0.33201 | -0.0386 + 0.31731 
-0.6826 + 0.45291 -0.0429 - 0.33201 -0.0386 - 0.31731 
-0.6826 - 0.45291 -0.1767 -0.1769 


Table 7. Full state feedback (partial distribution) with disturbance feedforward control law 
optimization results and performance 
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Figure 22. Simulation with full state partial distribution control and disturbance 
feedforward, sea state three (head seas) 


8. Integral Control 
This full state feedback with partial distribution was augmented with integral control 
on depth to remove the average depth error. As before, the integral control was done using 


the bow planes only. This results in the following control law: 


W — W commanded (1 04) 
D B К, | 0 0 M р у д M 
> ~ “commanded 
Os, О Кон № р 
< — commanded 
21 
|< ó (105) 


max 
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After a stable set of random gains was determined, the controller was optimized to 


minimize the deviation from the average depth. The formal optimization statement is as 





follows: 
Minimize: 
(106) 
F(Ky, Ku, Kis, Ky, Ks, H,F)— 
Subject to: 
real(eigenvalues(A.)) S Eqnay (107) 


This approach was used for each of the four sea state cases. For sea state three (head 
seas), the optimized response is shown іп Figure 23. The results of the four optimizations are 


shown in Table 8. 
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12.5543 











22.5268 
24.9900 






H/F (10? pounds) 20/0 20/0 20/0 20/0 
Mean Depth (feet) 55.05 55.31 55.21 
RMS Error (feet) 0.2106 0.533 0.4868 1.205 
LE 28 75 
0.0077 0 0.0003 0 


= 5 
26.7/4.6 21:2/-0.2 10.5/-2.6 


КЕ ' 
22.5268 
24.9900 
. 0 
0 












-0.0412 + 0.35871 
-0.0412 - 0.35871 
-0.1785 


-0.0412 + 
0.355 
-0.0412 - 0.35871 
-0.1785 
-0.0004 















15,22. ж Ре a) 
zi A ES 

EI eR UR. 
MO Мз. и; 2 
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% 
SN 
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-0.0412 + 0.35871 -0.0412 + 0.35871 
Mean Depth (feet) 7502 55.05 59:28 55.15 


-0.0412 - 0.35871 -0.0412 - 0.35871 


RMS Error (feet) 0.059 (28%) ОБО) 0.4274 (88%) 0.909 (75%) 
Maxtmum Error (feet) 0.248 (29%) 0.862 (50%) 1.16 (64%) 3.59 (105%) 


Eigenvalues -30.7579 
-4.6840 -1.8146 + 1.65551 -0.1526 + 0.32971 
-1.0562 -1.8146 - 1.65551 -0.1526 - 0.32971 
-0.1695 -0.6572 -0.1168 
-0.0005 0.00003 -0.0001 


Table 8. Full state feedback (partial distribution) integral control law optimization results and 
performance 
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Figure 23. Simulation with full state partial distribution feedback integral control, sea state 
three (head seas) 


D. FULL STATE FEEDBACK 
de Basic control 
The best control possible using state feedback should result from the use of all four 


states by each control. This results in the following control law: 


3b 


W — Wecommanded ( | 08) 


М | m s ] K, 2 К, 3 К, 4 Я = Q commanded 


m Ко Ку К, к Uu dn 


27 <commanded 


|< 5. (109) 


After a stable set of gains was determined using a linear quadratic regulator algorithm, 
the controller was optimized to minimize the deviation from the average depth. The formal 
optimization statement was: 


Minimize: 


(110) 





Subject to: 
real(eigenvalues(A.)) < E nax (111) 


This approach was used for each of the four sea state cases. For sea state three (head 


seas), the optimum response is shown in Figure 24. The results of the four optimizations are 


shown in Table 9. 
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H/F (10? pounds) 20.1/-0.9 18.0/-0.3 18.0/-0.33 19.9/0.0 


RMS Error (feet) 0.037 (40%) 0.2638 (73%) 0.2683 (76%) 1.24 (78%) 
Maximum Error (feet) 0.119 (45%) 0.961 (729) 0.961 (88%) 4.06 (89%) 


-1.0446 -1.0313 + 0.42181 
-0.2510 -0.3493 + 0.28811 | -0.3493 + 0.28811 -1.0313 - 0.42181 


-0.0617 + 0.52641 -0.3493 - 0.28811 -0.3493 - 0.28811 | -0.0936 + 0.07411 
-0.0617 - 0.52641 -0.2053 -0.2053 -0.0936 - 0.07411 





Table 9. Full state feedback control law optimization results and performance 
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Figure 24. Fullstate feedback optimized control simulation, sea state three (head seas) 


2% Disturbance feedforward 


The state feedback control law was modified to include disturbance feedforward. This 


results in the following control law: 


W — Wcommanded (1 1 2) 
д, [e Ку К\з a q — Q commanded "ext "rd F, 
QU K^ Ку Кэз Кол C= НЕЙ C Koa С, Koa Ма 
£ — £ commanded 
l| « à. (115) 
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After a stable set of gains was determined using a linear quadratic regulator algorithm, 
the controller was optimized to minimize the deviation from the average depth. The formal 


optimization statement Was: 


Minimize: 





A (114) 


2 
J (2 — £Z commanded ) dt 


0 


F(K u Кур, Куз, Куд) Ку, Кар, Код, Кра, Фе» Н, Р) = 


со > 
y 


Subject to: 


real(eigenvaluestA.)) S E max (115) 


This approach was used for each of the four sea state cases. For sea state three (head 
seas), the optimum response is shown in Figure 25. The results of the four optimizations are 


shown in Table 10. 


39 






Sea State / Direction 3/head 3/beam 
... ЗЕ: RARAS > ae Vs Tr A AGAR 
Initial Values АСА З AE GAS у EP 37,7070, AS 
EII SAS AS Fes Sas oa Shel lC pa 
OA қр МҰҚАН Еле e iouis oy НО 
"ОРИК ЧИЛ ГЛ ЛА ИЕ RSI ди, р" о ЗБ рву placer sue pituy E= Аче у 
KT 


Eigenvalues -0.8859 -0.8859 -0.8859 -0.8859 
-0.2854 -0.2854 -0.2854 -0.2854 

-0.1630 + 0.22471 -0.1630 + 0.22471 | -0.1630 + 0.22471 | -0.1630 + 0.22471 

-0.1630 - 0.22471 -0.1630 - 0.22471 -0.1630 - 0.22471 -0.1630 - 0.22471 


š У 5 
25, Mag 2 ⁄ + 2: сер Am Y y T к” 
АЕ e e V ME E posee A; T Ee ddl f ^ Р 
ФРА кез бе. [рр чол P AE Co Ter "Line ES ОН ROS к RE Ўз st i 
Қайғысы роза Z PAD SEES! 34 ы A A e E Я 6 AA 
K : | : s | | -3.94 
: е | | 600.1 
| ; | ; 27.66 48.6 | 449.5 


0., (rad/ sec) 1.4046 0.983 1.739 
H/F (10? pounds) 13.5/-0.8 19.16/-0.2134 19.0/-0.1 18.8/-0.1 






RMS Error (feet) 0.0928 (65%) 0.4121 (31%) 0.400 (36%) 0.792 (33%) 


Maximum Error (feet) 0.2852 (58%) 1.62 (47%) 1.117 (101%) 2.23 (41%) 


-/ 0591 -0.8095 + 0.52681 -1.1927 -0.0898 + 0.92971 
-4.6280 -0.8095 - 0.52681 -0.2499 + 0.33121 -0.0898 - 0.92971 


:0005225::0:15291 -1.1529 -0.2499 - 0.33121 -1.1834 


2001522 - 05193291 -0.0315 -0.0618 -0.6517 





Table 10. Full state feedback control law with disturbance feedforward optimization results 
and performance 
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Figure 25. Full state feedback control with disturbance feedforward optimized control 
simulation, sea state three (head seas) 


B. Integral control 

This full state feedback with partial distribution was augmented with integral control 
on depth to remove the average depth error. Since the bow planes are principally used for the 
control of depth, the integral control was done using the bow planes only. This results in the 


following control law: 


W — W commanded (1 1 6) 
4 — 4 commanded 
Op | | Ki, Ко Ку Ка К; 0 5 EE 
5 = K K K K K Т “ commanded 
sp 21 22 23 24 25 
€ — Ccommanded 
L < 
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After a stable set of random gains was determined, the controller was optimized to 


minimize the deviation from the average depth. The formal optimization statement is: 


Minimize: 





, (118) 


2 
| (2- C commanded ) dt 


0 
ЕСК, Ку, Куз, Куа, Ку5, К, Ко, К, К, Къ5. Н, Е) = 


ry 


Subject to: 


real(eigenvalues(A.)) S Emax (119) 


This approach was used for each of the four sea state cases. For sea state three (head 
seas), the optimum response is shown in Figure 26. The results of the four optimizations are 


shown in Table 11. 
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-0.1652 + 0.21531 
-0.1652 - 0.21531 
-0.122 


256.28 
84.99 
0.1753 


-0.0128 


20 759 


-17.3244 
-0.7017 
-0.1942 + 0.42961 


-0.1942 - 0.42961 
-0.0076 


24.2736 





-0.1652 + 0.21531 
-0.1652 - 0.21531 
-0.122 


трн.” Мо ИАА, ОТ 
Е 
^ я y ot UA ЛЕЙ £ Y 
SPE ae ауы 


-0.1652 + 021531 
-0.1652 - 0.21531 


75561 |715: 


153.309 
39.280 
-0.0550 
0.0069 


Mean Depth (feet) 55.00 


-0.8739 
-0.1541 + 0.14411 
-0.1541 - 0.14411 

-0.2567 

-0.0375 


-0.2086 + 0.16281 


-0.2086 - 0.16281 
-0.2056 
-0.0208 





4/beam 


UN Y з Багд рари ^ 
Pe ONO ata tesa 


-0.1652 + 0.21531 
-0.1652 - 0.21531 


48.260 
-0.0327 
-0.0105 





19.9/0.0 
1.96 (79%) 


6.88 (99%) 
-1.0690 
-0.1214 + 0.33271 
-0.1214 - 0.3327 
-0.2428 
-0.0012 





Table 11. Full state feedback integral control law optimization results and performance 
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Figure 26. Optimized full state feedback with integral control simulation, sea state three 
(head seas) 


E. CONCLUDING REMARKS 

For each case of feedback control, the degree of control achieved generally improved 
by the additional state feedbacks. Full distribution of each state to both controls further 
reduced the error. Table 12 provides a summary of the optimizations performed, and the RMS 
error of each one. 

Changes in the optimized trim and control law in all cases varied substantially with 
changes in sea state or direction. This is consistent with operational experience. 

Each controller was optimized with only the goal of minimizing the mean square of 
the depth error. This resulted in large gains and excessive control effort. In addition, large 
rates of control were experienced. This would be detrimental for actual submarine operations, 


as there are rate limits associated with the control surfaces. These limits come from the sizing 


64 


of the hydraulic plants which drive the planes, and operation concerns related to plane induced 
noise. 

Some improvements in depthkeeping were achieved by the feedforward of the 
disturbance forces. This is in spite of the feedforward being based on a steady state response 


to a constant disturbance. 









Sea State/Direction 3/Head 
Control Scheme Е OR 










Table 12. Optimized RMS error (feet) of state feedback control schemes 
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V. SLIDING MODE CONTROL 


A. INTRODUCTION 
Í. Overview of MIMO sliding mode control 
The controller design starts with a standard linear state space representation: 


х= Ах+ Ви (120) 


where: 


A € R”™™ | state matrix 
B eV". control matrix 
x € U" state vector 


ue R™! control vector 


The sliding mode control law, 4, is composed of two main parts: 


u=u+u (121) 
The first part, & , 1s a linear feedback based on the linear representation given by Equation 
(120). The second part, @ , are nonlinear feedbacks with their signs switching depending on 
the relationship of the system states to the sliding surfaces. The sliding surfaces are hyper 
planes in the state space, with one for each control. They are defined by: 


п) Ет (122) 


Where: 


SER” 


To determine the nonlinear feedback functions, the concept of Liapunov stability is 
used. The Liapunov function is taken as: 


(123) 


V(x) = 5 (7 + 02 *..*0;) 
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Asymptotic system stability is guaranteed provided that V(x) is a positive definite function, 
that 15: 


V()=0,0, +0,0,+..+0,0, <0 (124) 


Equation (124) is satisfied if: 
0,0; «0 for 7 71 to s (123 


Equation (125) can be rewritten as: 


O; 7 —rsign(o;) (126) 


where T, is a positive gain parameter for the 7^ sliding surface. Equation (126) can be rewritten 


after substututing the tume derivatuve of Equation (122) and Equation (120) as: 


7 sign(o,) (127) 
Np Sign(O> ) 
51 (Ах + Ви) = – 
nN, Sign(O,, | 
Solving Equation (127) for z yields: 
ї5ї&п(с\) (128) 
115 sign(O») 


u == (sip s n = c py 


1, sign(o, ) 


Equation (128) can be rewritten in matrix form as: 


u — Kx K,nsign(o) (129) 


where: 
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Kes s py СТО 


Koi Ep 


Equation (129) is identical in form to Equation (121), with a linear state feedback and a 
nonlinear switching term. 

For the decomposition in Equation (128), it is required that the closed loop stability 
matrix, (4-BK), have n zero eigenvalues. The sliding surfaces are the left eigenvectors resulting 


from the zero eigenvalues. 


D. Utkin's method for MIMO sliding mode control law design 
Determination of the sliding surfaces can be difficult, especially with a MIMO system. 
One technique for this is proposed by Utkin (1977). For this technique to be applied, the B 


matrix of the state space system (Equation (120)) must be of the form: 


р ІН (130) 


where: 


B, e 9i" 


0 € R (т-п)хп 


For the MIMO cases of vertical plane depth control used in this thesis, this was the case. Por 
the stern planes only control examples, a QR factorization was applied, to transform the state 
space system into this form. 

Given that the B matrix is of the form defined in Equation (130), Equation (120) can 


be decomposed into the following: 


И Арх СВ (1995 


Ху = Ах) + А;2Х2 (132) 
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where: 
хе" 
ОА 
Е 
Е 
ЕО 


(m-n)x(m-n) 
Az, ER 


The sliding surfaces become: 


SAS 0 (158) 


where: 
5! ev" 


55) c gg x(m-n) 


Because the sliding surfaces are the left eigenvectors of the z zero eigenvalues, $ can 


be set to the identity matrix without loss of generality. Substitution of Equation (133) into 
Equation (148) leads to: 


С = X) + Š; x = —ñsign(xi + S; x>) (134) 


Thsign(o)) (135) 
15 sign(O5) 


uz "n en + SZ A, )x) + (A) TAT | n 


| n,sign(o, ) 
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When the system 1 on the sliding surfaces, Equation (133) can be used to solve for x, 


in terms of x,, resulting in: 


е у 2 (130) 
Substitution of this result into Equation (132) results in: 
x, AAT SEE (137) 


Equation (137) 1s the set of independent equations that the non zero eigenvalues for 
the control result from. For the application of pole placement algorithms to determine the 


sliding surfaces, it only has to be recognized that it is in the standard state space format: 


у = Ау+ ВКу (138) 


with: 


Once $57 is determined, the control law can be determined by substituting it into Equation 


(135). 


Utkin's technique also allows for the use of linear quadratic regulator methods for 


determination of the sliding surfaces. (Utkin 1977) For this method, it is desired to minimize a 


quadratic performance index: 


= 139 
| = 2 [xt Qxdi 0 
- 0 


pn 


where Q is a positive definite weighting matrix. By partitioning O in the same manner as A 


was partitioned for equations (131)and (132), Equation (139) can be rewritten as: 


Ef rur T T T p 
127] 6 Qux tX rx ODE Ош 
0 
ОК 
17 ! (141) 
І =– 070 xy tv" Ov" at 
0 
шеге: 
Q - 05 - 01010» 
А = As An Qu Qro 
v = x *Qn OQ 
With the system on the sliding surfaces, Equation (140) can be rewritten as: 
А = Hw T Av (142) 
Equations (141) and (142) are in a recognizable form for the application of any 
convenient linear quadratic regulator solution. The Hamiltonian is: 
(143) 


* | ж 
H = р'(А х, + Аду) > (2 Q Оу) 
The algebraic Riccati equation 1s: 


We 


(А) КЕХА —КАОГ АКО =0 (144) 


The solution of Equation (144) results in: 


v2-QuAlku (145) 


This result is used with the definition for v from Equation (141) to provide the relationship 


between x and x». This results in the sliding surface: 


52 = Оп (Ор + Ад) (140) 


With 57 determined, the control law can be determined by substitution into Equation (135). 


This sliding mode LQR controller design was implemented in a MATLAB? function 
SMLOR.M which included provisions for a OR factorization for the cases when the B matrix 
was not of the form given by Equation (130). SMLOR.M is included in Appendix A. 

3. Control of chatter 

One undesirable aspect of sliding mode control is the chatter induced by the nonlinear 
switching term near the sliding surfaces. The nonlinear switching term in the sliding mode 
control can cause control chatter when the system 1s near a sliding surface. One way to reduce 
the chatter is to use the concept of a boundary layer around the sliding surface. First, a satszgn 


function is defined: 


Lx»1 (147) 
satsigentx)=3x -1<x<l 


—-]lx«-1 


A boundary layer of thickness ¢ around each sliding surface, is applied using the satsign 


function: 


> 


| [5 | (148) 

С; = satsign| —— 
0; 

Equation (147) can be used to replace the szgv function used in Equation (129) with no change 

in the asymptotic stability of the system. There are some effects, however, because the 

dynamics near the sliding surface are not the same as the closed loop dynamics which exist 


when the system is on the sliding surface. The final control law is: 


149 
u= Kx+ клава Z e 


Block info and requirements 


Saturation Sumi PD Control Remove state vector 


commanded 


Feed Forward 





Figure 27. SIMULINK® model sliding mode controller 


Equation (149) was implemented as a SIMULINK? model, shown in Figure 27. 


B. SIMO SLIDING MODE CONTROL RESPONSE TO DISTURBANCES 
When applied to vertical plane submarine control, sliding mode control has several 

nuances which are not obvious from inspection of the governing equations. In order to 

illustrate these, the performance of sliding mode control will be explored through several 


example cases. To keep the analytic derivations simple, the cases worked will be done with 
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stern planes control only. The general concepts, however, will be applicable to stern and bow 
planes control. 

The response of sliding mode control to force and moment disturbances 1 
fundamental to its application to submarine control in the vertical plane. These force and 
moment disturbances could result from a variety of situations. Examples of these include out 
of trim conditions, free surface effects, and wave forces. 

]. Basic sliding mode disturbance response 

The first case will use a basic sliding mode control law, that is one without a 


feedforward term or integral control. The resulting control law is: 


O=S w+S,q+S,0+S, 2 (150) 
O (151) 

ô = Kw + K,q + КО + о) 
6| « à... (152) 


a) Linear analysis steady state 

Assuming that the sliding mode control is not saturated and that the control 
deflection 1s less than the maximum, the submarine should reach an equilibrium state under the 
action of steady force and or moment disturbances. The linear equations of motion in the 


vertical plane with one control are: 


w 5 ayuw- a,4uq t a,30 bu, Ó +F; (153) 
¿=a,uw+a,4q +00 +b,u ô+ M, (154) 
йер (155) 

i-w-u8 (156) 

xv (157) 


Equations (150) through (157) can be solved for the following steady state condition: 


Э 


Pub (158) 


w = 
2 2 
b,a,u а ban DIE 


ХА 


о (159) 


_ My) (160) 
= b,a u? +b,a,-ba,u? -ba 


a, + K, b, +4, фи К 0b,S, +u*nK Qb,S, + K u*b, (161) 
u’ K_S,0(b,a,,u? +b,a,,—b a, u* —b,a,,) 

Куи*Ь, + ауи* + ау, + и?тК ,фЬ,5| + Kub, +u nK, ob 5, 

u*K,S,p(b,a, u? +b,a,,—b,a,u* —b,a,,) 


7 = Г d 


ce 


T (a, u° +a,)+ M, (баци: — a) (162) 


2 2 2 
u“ (b,a, u +b,a,, —b,a, u” —b,a,3) 


NS 


If the sliding mode is just saturated, that is |o / 9| 2 1, the system will still be 
stable, however the gain parametern will be at a critical value. If further reduced, the system 
will be unstable. Assuming the sliding mode is just saturated this critical value, 7],,, , can be 


determined. 


(163) 


али’ +a, + Kyu) b, * K u*b, 
d 
K u*(b,a, u* +b,a,,—b,a, u* —b,a,,) 
2 2 3 
йш а кк Km 


d 2 2 2 


T = 


For cases with the gain parameter 7] less than critical, a steady state Z results, 
along with steady state values in all states other than z. By assuming a steady state 2 апа 


solving Equations (150) through (157) for equilibrium, the following equations result: 


Ñ и“ Куба —-bjay) * Ку (азі — 4230) ) + 413421 — 423411 
(165) 


3 
| Je, рін Мең Өн AOS 
EE A AA c 
а u(b uK; + a;,) 
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Wyr = UO d, (166) 


9 = 0 (167) 
5 = (4, Uw,, + 2130 + F) (168) 
AN bu 


Inspection of Equation (164) may cause the reader to incorrectly assume that a 
nonzero steady state value of @ will exist in the absence of disturbances. This is not the case, 
however, because the use of Equation (164) implies the lack of a steady state z , and therefore 
a disturbance. 

b) Nonlinear analysis steady state 

An analysis similar to that conducted on the linear equations can be conducted 
to determine the system steady state response under a constant disturbance. The nonlinear 


equations of motion in the vertical plane with one control are: 


w- ay uw a4uq * aj, sin(0) +b,u28 + F, cos(0) + e,,9* +e,,qw (169) 
д = ауну +аууид + ауу 5\п(Ө) - bu* 6 - M, cos(0) - e,q^ ey, qw (170) 
д=а ШЕШІ) 

¿ = wcos(0) — usin(0) (122) 

х = wsin(0) * ucos(0) (173) 


Once again, the basic sliding mode control of Equations (150) and (151) is 
used. Assuming that steady state 1s possible (sliding mode control 1s not saturated, and the 
control deflection is less than the maximum) Equations (169) through (173) subject to the basic 


sliding mode control law can be reduced to the following nonlinear equation: 


aj, ѕіп(Ө)и?Ь, + ауз ѕіп(Ө)соѕ(0)Ь, + Е, соѕ?(0)Ь, – (174) 
ba, sin(9)u? — ba; sin(0)cos(0) - bj M, cos^(0) 2 O 


Tm 


Solution of Equation (174) yields the steady state value of 0. It is of note that 
Equation (174) is independent of the control law. The correct root is readily determined by 
using the value closest to the linear analysis (Equation (160)). Substitution of this value into 


Equations (169) through (173) yields following steady state results: 


w,, =utan(@, ) (175) 
q. = 0 (176) 
5 = (Quuw, +apsin(0,)+ F,) (177) 
NN bu? 
-E (8, - Ko, - K0,) - Sov, 50, (178) 
С = 1 - 
Я 2 


The steady state value of z given by Equation (178) Is dependent upon the 
control law gains. Following the example of the linear analysis, a critical value of the gain 


parameter, т} can be determined. 


бу – Куру – Куда (172) 


Merit = x K 





5 


Unlike the linear analysis, the critical value of the gain parameter is not just a 
linear combination of the disturbance forces, but rather requires nonlinear solution for each 
possible case. 


If |o / $| » 1 the sliding mode control will be saturated, and a nonzero steady 


state Z.. will exist. For this case, the control law control law reduces to: 


ó 2 Kw Ка + КӨ + nK,sign(co) (180) 


Given the steady state final condition in all state variables with the exception of 
z, Equations (169) through (173) and (180) can be reduced to finding a root of the following 


nonlinear equation: 
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u* (bj K4,a5,0 * b, K,a5, - b Kyaj0 - b; K5a,,) + usin(@)(a,3b) K, — aj4b Kj) p (181) 
byuK, +4» 

F, cos(0)Mb7uK, + a7,)- M , cos(8 Mb uK, +a,,) +sin(9)(a¡34,, — 493011) — 
b,uK, +a, l 


0 


This can be accomplished using the using the linear value of 0,, from Equation (164) as the 


initial guess. It follows that: 


m cos? (0,.) + a); sin(@,.)cos(@,, ) (182) 
fsx 7 u(b,uK, + a,,) 
u*(a,, sin(@,,) +b,uK, sin(@,,) +b, K30,, cos(@,,) - b, Kr] cos(0,,)) 
E u(byuK; +a,,) 
Wss E tan(6,, ) Us 25 (1 83) 
ды E К, Wss + КзӨ,, T nK,sign(o) (1 94) 
с) Disturbance response simulation with basic sliding mode 


These results can then be applied to the SUBOFF hydrodynamic coefficients. 
For the modified coefficients used in this thesis, the linear state space system at six knots with 


stern planes only is: 


У] [-0.0179 37101 00196 0lw] [-01009 F, (185) 
ġ | | 00006 -00680 -00034 O|q| | -.0027 M, 
0 | о l 0 ole wee | 
E | 0 101269 0(: 0 0 


A sliding mode control law is determined using Utkin's method. After some 
experimentation, the diagonal of the minimization matrix Q was selected as Qu = 100, Qz = 


100, Q33 = 100, Qu = 1. This yielded the following control law: 


с = 1.0и/- 6755924 —15.40820 + 0.08305 (186) 


Do 


б (462) 
Ô = -0.2803w + 84.8572q +7.06120 + K nsatsign - 


[6] < 04 (188) 


A moment disturbance was chosen to be controllable but give a nontrivial 
response. The application of a force of five thousand pounds and moment of 4,573 thousand 


foot pounds resulted in a pure angular acceleration , M, =0.001 radians/second?. For this state 
space system, control law, and disturbance set, the value of 7,,;, was (0.0661). Application of 


the nonlinear solutions yielded the steady state solutions in Table 13. 


о | = om | wes wes [ ese 
[o oum [oy oem pes | sma sem 
Lu poe pompe ems - 
[3 o9 | 9 | таи | пат | за | 3 — 
| аза ov pes ues | wm | 3 — 


Table 13. Steady state nonlinear solutions for M, 20.001 radians/second? 









The system transient response was simulated using the RK45 function of the 
SIMULINK® toolbox. Figure 28 shows the resulting paths. The response 1s given for six 


values of the ratio of n and N. As expected, for values of n less ап спіса! е сопко] 


law is unable to maintain a steady depth. 
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Figure 28. Nonlinear simulation of vertical plane response to a pure moment disturbance 


The calculations and simulations were repeated for a pure force disturbance. A 
force of 43 thousand pounds and a moment of 220.4 foot thousand pounds resulted in a pure 
vertical acceleration of 0.005 feet/second?. For this state space system, control law, and 


disturbance set, the value of 7,,,, was (0.0466). 


w (Feet / Sec) q, (Rad / sec)| @,,( Degrees) O, (Radians) | 2, (Feet! Sec 


1.5724 5.4888 0.2357 0.5966 






AAA s — 

Lo ew D$ | mer ras 
[ow pue | s se p mam sme | o 
[3 ue | s wee mer oor 0 
[we p» wee was [wm | 3 — 


Table 14. Steady state nonlinear solutions [ог Е, = 0.05 
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Figure 29. Nonlinear simulation of vertical plane response to a pure force disturbance 


d) Disturbance response with sliding mode feedforward control 

To eliminate the steady state error induced by a disturbance, several techniques 
may be employed. Using a feedforward function in the sliding mode control law is one 
technique. A feedforward term works by using knowledge of the external disturbance and 
using applying some degree of control effort. This provides a steady state control effort to 
oppose the disturbance without a steady state error. This approach ts limited as it requires one 
control per zero error state and may be limited in other ways. A feedforward term can be 


added to the sliding mode control law by changing the sliding surface to the following: 


OG =S,w+S,q+S,6+S, 2+5S, (189) 


The value of S, is such that it will equal the control effort that is applied by the 


steady state quantity that is desired to be zeroed. Since the main priority is to obtain zero 
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depth error, S$, is selected to equal the control effort introduced by the previously calculated 


value of steady state depth error: 


S = 8 = K WS К.О, )- SW. = 5-6 (1 90) 


ASN 


A 


The steady state values needed for Equation (190) can be determined from 
linear or nonlinear analysis, although the linear analysis will result in a non-zero depth error. 


The linear analysis gives the following: 


$ =СМ,+СЕ, (191) 


where: 





QS, um (192) 


- ? e + e Kb, + Kub | 
WIxNCLSU u | 54 54 


on 2 2 


193 
С SU Ks, Кв, Зи ( ) 


DUNT 


2 2 
аи b, + аЬ, -bau — һа» 


5 ES 
С, = 4 Л] 


83 


Depth, Feet 





0 500 1000 1500 2000 2500 
X axis, Feet 


Figure 30. Nonlinear simulation of vertical plane response to a pure moment 
disturbance with a feedforward term based on nonlinear steady state 
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Figure 31. Nonlinear simulation of vertical plane response to a pure force disturbance with 
a feedforward term based on nonlinear steady state 
e) Disturbance response with sliding mode integral control 
Another means of eliminating steady state error is by the use of an integral 
control term. To accomplish this, an additional equation is added to the state space 


e 


representation. 


opes (194) 


This forces a zero steady state error in z , although there are some additional 
considerations with the use of integral control. The resulting control law, with the additional 
State, 15: 


mc (196) 
О = Ант АА О KZ THES CHE MM 
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|6] < 0.4 (197) 


Based on inspection of Equations (195) and (196) several conclusions can be 
drawn. First, because z is included in the proportional portion of the control law, values of the 
gain parameter which are less than critical will result in a steady state error in z for any 
controllable disturbance. Second, if there is a steady state error in z , the magnitude of the 


integral term, z, , will tend to infinity. This can cause problems with changing conditions or 


pathkeeping as it will delay the control response to other errors. 

At a condition of steady state, with the gain parameter greater than critical, the 
steady state error in zis zero. Because of this, the previously calculated values for steady state 
pitch angle, heave, and control deflection are still valid (Equations (174) through (177)). 
Moreover, because any controllable disturbance will result in a steady state condition these 
values apply for cases where the gain parameter is less than critical. 

For conditions where the gain parameter is less than critical, the resulting 
control law 15: 


ó = Кү + Куа + K30+ K¿2+NK,sign(0) (198) 


And the steady state value of z 1s: 


ô.. — Kiwa — Коб, — K,sign(O) (ДЕР) 


55 = K 
4 


Expressed in linear state space form for control law design, the representation 


for the SUBOFF at six knots 1s: 


и] [0.0179 37101 00196 O0 Ol|[w] [-0.1009 n (200) 
d 00006 -0.068 -0.0034 0 Ol q | |-0.0027 M, 
0|l-| O 1 0 О | о 1a 41) 0 
i 1 0 -—-101269 0 Ol = 0 0 
AME Ee 0 0 | Olz, 0 0 
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Utkin's method was applied to obtain a suitable controllaw. After some 
experimentation, the diagonal of the minimization matrix Q was selected as Qui = 100, Qz = 


100, Озз = 100, Qu = 1, Qss — 0.01. This yielded the following control law: 


с = 1.0% – 49.81584 – 1277250 + 0.1050: + 0.00352, (201) 
5 (202) 

д = —1.6009w + 161.0589q + 24.81250 — 0.099 Iz + K .nsatsign n 
5|< 04 (203) 


Simulations of the SUBOFF under sliding mode integral control with а 
moment disturbance are given in Figure 32. For this state space system, control law, and 


disturbance set the value of 7,,;, was (0.0484). Shown are five different values the ratio of the 


gain parameter to the critical gain parameter. For values of the ratio larger than about two, the 
system exhibited excessive oscillation before settling to zero depth error. 

simulations of the SUBOFF under sliding mode integral control with a force 
disturbance are given in Figure 33. For this state space system, control law, and disturbance set 


the value of 7,,;, was (0.0466). Shown are six different values the ratio of the gain parameter to 


the critical gain parameter. For values of the ratio larger than about two, the system exhibited 


excessive oscillation before settling to zero depth error. 


1) Sliding mode disturbance response conclusions 

Submarine vertical plane depth control using sliding mode control can be 
effectively achieved in the presence of disturbances. Sliding mode control is similar to linear 
state feedback in that the an external disturbance will result in a steady state depth error. 
However, if the gain parameters of the sliding mode control are not properly selected, a loss of 
depth control can occur. 

Steady state error can be dealt with using feedforward or integral control. 
Integral control has several advantages. Application of integral control does not require 
knowledge of the disturbance. However, if the gain parameter is too small, windup of the 


integral term occurs. If the gain parameter is too large, excessive oscillations can be 
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introduced. The greatest advantage of integral control demonstrated was that gain parameters 


less than critical did not result in loss of depth stability. 
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Figure 32. Nonlinear simulation of moment disturbance using sliding mode integral 
control 


Feedforward control exhibited good disturbance compensation, assuming that 
the disturbances were measurable. Given the disturbances, feedforward values can be 
determined based upon the nonlinear equations of motion, requiring periodic nonlinear root 
finding, or upon the linear solution. Because of the computational expense of obtaining the 
nonlinear solution and the expected error in the hydrodynamic coefficients, a linear steady state 


solution 15 appropriate for feedforward computation. 
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Figure 33. Nonlinear simulation of force disturbance using sliding mode integral control 


C. MIMO SLIDING MODE CONTROL AT PERISCOPE DEPTH 

|; Introduction 

The purpose of using sliding mode control was to provide an alternate means of 
control which relied upon all the system states. The robust characteristics of sliding mode 
control were thought to be a good approximation to the human operators. 

At periscope depth, experience dictates several desired conditions. First, the ship is 
trimmed heavy to counter the steady wave forces. Even more weight 1s brought on after this 
point to allow for a constant small positive trim angle, of several degrees. This provides 
reserve ballast which is made available by reducing the trim angle. Finally, in sea state three, it 


should be very possible to maintain depth within one foot of ordered depth. 


89 


The equations of motion used for this section are the nonlinear equations of submarine 
motion in the vertical plane. They are different from the equations used previously in this 
chapter, as they include both bow and stern planes. Also the force and moment disturbances 
used represent not only constant disturbances, like ships trim, but time varying wave forces as 


well. Repeated for convenience, the equations are: 


у = аууиу t a,4uq ta, sin(0) 4 b, E tb ,u*6, * F, cos(0) +e, ‘aa +€,,qw (204) 
à 2 ag uw t a55uq * a3, sin(0) 4- b, u26, Бш б ЕМ Cos(0) + 2,14 теп (205) 
ó-q (206 
z 2 wcos(0) — usin(0) (201 
x 2 wsin(0) + ucos(@) (208) 
2. Basic sliding mode controller 


The sliding mode controller 1s of the same form as before, although with the 
introduction of MIMO control some of the scalar terms become vectors or matrices. 'The 


form of the basic sliding mode controller 15: 


O^ — S, W + 5224 + 5-30 F S54 2 (210) 

mic mics (211) 
Ôp = Kuw + К124 + К\30 +11 К, satsign Lx * 15; K,,,satsign ГЕ 
1 2 

NE NE (22) 
б, = K,¡w+ Kng + Куд + К, satsign E + 2K,,,Satsign 7" 
l 2 

б, | S О тах (213) 

5, ES О (214) 
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As an initial attempt, Utkin's method was used to determine a sliding mode control 
law. After some experimentation, the diagonal of the minimization matrix Q was selected as 


Оп = 100, 022 = 100, Озз = 100, О = 1. This resulted in the following control matrices: 


[396585 -63281 -40454 0 (215) 
 |-246916 42965 25027 0 
l 0 (216) 
1 0 ] 
|—-0085 17281 


0.0170 -00985 


y 


10.9711 —405.06 (217) 
3083 25210 


The values of the gain and boundary layer thickness parameters were taken as: 
"ed (218) 


ф = ф, = 1 (210) 


A SIMULINK® model was developed to incorporate the submarine dynamics of 
Chapter Il, the wave forces of Chapter III, and the MIMO sliding mode control law. Also 
included was the trim model from Chapter IV. 

The model was used to simulate a step change in commanded depth from 140 feet to 
50 feet in depth. To provide some realism in the trim condition, the submarine was trimmed 
to 25 thousand pounds heavy, with no moment correction. Wave force values for sea state 
three were used, with a relative heading of 180 degrees (head seas). Figure 35 shows the 


resulting path taken by the submarine. 
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Figure 34. SIMULINK® model of submarine with wave forces and trim 


At first glance, this control scheme fulfills most of the desired characteristics of a 
submarine depth controller. Periscope depth was achieved with no overshoot, and reasonable 
depth control was maintained. The trim condition was selected so that a steady state positive 
trim angle would exist at periscope depth. During the depth change, the maximum trim angle 
achieved was about ten degrees, which is also very consistent with actual submarine practice. 

Inspection of Figure 36 shows some problems with this particular controller. The 
application of control effort was excessive. The main reason for this was the high frequency 


variations in wand q induced by the wave forces. Because of the combination of wave forces 


and trim, the commanded depth was not achieved. The average depth at periscope depth was 


50.75 feet, as opposed to the commanded depth of 50 feet. 
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Figure 35. Basic sliding mode performance, step change approach to PD 


Although the performance for the sea state three, head seas was adequate with this 
controller, it did not perform well with the other sea states or headings. Because of this, it was 
decided to use this control as a starting point for a performance optimization for each of the 


four sea state and direction cases available. 


As was done for state feedback control, the MATLAB® CONSTR function was used 
to perform the optimizations. To provide a general set of design variables, the sliding mode 
linear quadratic regulator program was not used to determine the control law at each step of 


the optimization. Instead, the sliding surface itself was varied to change the control law. 
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Figure 36. State parameters for basic sliding mode approach to periscope depth 


The formal optimization statement 15: 
Minimize: 
(220) 





where: 





Z mean E 


H = Ballast added to center of buoyancy, thousands of pounds 
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F — Ballast shifted from forward to aft, thousands of pounds 


Subject to: 


real(eigenvalues( A5; — AOS E (221) 


Deviation from the mean value of depth vice the commanded was used because of the 


expected average depth error. 


This approach was used for each of the four sea state cases. For sea state three (head 
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Figure 37. Simulation with basic sliding mode control in sea state three (head sea direction) 
seas), the optimized response is shown in Figure 37. The results of the four optimizations are 


shown in Table 15. For the RMS error and maximum error, the optimized values are given, 
along with their percentage of the initial values. In all cases, use of the optimization resulted in 
reduction of the mean square depth error (measured from the average depth). Reduction of 


the maximum error was also achieved. 
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Sea State / Direction 
Initial Values e 
S 


0 0 0 0 














0 0 0 0 


Sliding surface -0.8725 + 0.50631 -0.8725 + 0.50631 -0.6981 + 0.48201 -0.8726 + 0.41721 
eigenvalues -0.8725 - 0.50631 -0.8725 - 0.50631 -0.6981 - 0.48201 -0.8726 - 0.41721 


жү У, va о z s he 
Optimized Values A rac 
КО АРА МИ x 
RARAS 
EA ка š $t 
VU VC ER ФЕР АУА,» ENS YE 4 


Е 
0 0 0 
TA Е . : 

























0 0 0 
0.77 (26%) 2.36 (29%) 1.92 (46%) 


Sliding surface -0.8725 + 0.6945 -0.6968 + 0.54341 -0.8726 + 0.41731 





eigenvalues -08725 06945 -0.6968 - 0.54341 -0.8726 - 0.41731 


Table 15. Optimized basic sliding mode control law results and performance 
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ол Disturbance feedforward 


The sliding mode control can be implemented with a disturbance feedforward to 


correct average depth error. This can be implemented inside or outside of the sliding surface. 


For this example, the disturbance feedforward was implemented inside the sliding surface 


calculation. 


Оо 


(222) 


Assuming that neither sliding surface 1s saturated, that control deflection 1s within 


limits, and using a linear analysis, the steady state value of the depth error can be written as a 


linear combination of the force and moment disturbances (Appendix B). 


Lys 7 *commanded 7 Ci A E С, Ма 


(223) 


To eliminate the depth error, it is necessary to apply the same amount of control effort that the 


steady state error provides within the sliding surface: 


541 
5 p= (or — É commanded ) A 


This results in the following control law: 


Күз 
Кз 


0 
0 


| 


n W commanded 


Q — Q commanded 
0-0 


£ — Z commanded 


commanded 


W — W commanded 
q — q commanded 
0 > mania 


< —  commanded 


ó| < ó... 


n К. 1 
К, 


E 


+ 
С15 4 


where the force and moment disturbances are filtered. 


Оў 


(224) 
| (2. (225) 
7), satsign| — 
s Qi 
| | З >) 
1]; satsign 5 
(226) 
ЖИ 
C534; | Ма 
(227) 


The initial sliding surface and gains from the basic sliding mode control law was used 


as the starting point for optimization. The formal optimization statement is: 





Т (228) 


B 2 
| (5 — C commanded ) dt 


0 
F(5311,535,54,,545,0,,, 1, 15, H, F) — [ 
] 


Subject to: 


real(eigenvalues( Aj; — hoo e (229) 


This approach was used for each of the four sea state cases. For sea state three (head 
seas), the optimized response is shown in Figure 38. The results of the four optimizations are 
shown in Table 16. For the RMS error and maximum error, the optimized values are given, 


along with their percentage of the initial values. 
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Sea State/Direction 3/head 4/head 
e Muro 


Initial Values % = --- 
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0 0 0 0 


0.05/0.05 0.05 /0.05 


0.05/0.05 0.05/0.05 


Sliding surface -0.8725 + 0.50631 -0.8725 + 0.50631 -0.8725 + 0.50631 -0.8725 + 0.50631 
-0.8725 - 0.50631 -0.8725 - 0.50631 -0.8725 - 0.50631 -0.8725 - 0.50631 


x^ ien] c 
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-0.0981 1.7281 -0.1255 

0.0166 0.096 0.0400 | -0.0473 -0.0240 
-24.08 -11.70 -5.992 
429.65 530.65 467.3 
244.1 118.74 60.9 


0.05/0.05 0.0542/0.0578 0.1853/0.0544 0.0596/0.0548 
ЕЕЕ (105 pounds) 20/0 20:27:04 20.7 /1.5 20.1 /0.0 
RMS Error (feet) 0.405 (97%) 1722 (81%) 1.01 (1 00%) 1.99 (28%) 

Maximum Error (feet) 1.36 (111%) 3.20 (83%) 2.52 (1270) 4.93 (26%) 


Sliding surface -0.8723 + 0.48131 
-0.8723 - 0.48131 


0 0 0 


0 0 0 0 0 
















eigenvalues 





Table 16. Optimized sliding mode control with disturbance feedforward results and 
performance 
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Figure 38. Simulation using sliding mode control with disturbance feedforward 
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4. Integral control 
The sliding mode control law can be augmented with integral control on depth to 


remove the average depth error. This results in the following control law: 


W — Weommanded (2 30) 
С 
4- Q commanded satsign| —- 
HEN Ko Ki Ка I M К, Ку, L $ Es 
жеу ~ “commanded 
D Ко К > Кз Кд 0 К К. . | 02 | 
€ — Ccommanded T» satsign| ^. 
Ó 
2] 
] 0 | W — Wcommanded (25 1) 
0 ] 774 соттапае 
с, 4—9, ded 
= 34i 312 O = РУ, 
62 
Sa] 542 < — Ccommanded 
[Ss 952 | | ©] 
6] < 5 (232) 


тах 


After a stable set of gains was determined, the controller was optimized to minimize 


the deviation from the commanded depth. The formal optimization statement is: 










Minimize: 
p (239) 
| (Z= commanded I 
F(83,,835,541,545,555,355 Th T5, H, FP) — : г, 
Subject to: 
real(eigenvalues( Aj; — DE a (234) 


This approach was used for each of the four sea state cases. For sea state three (head 
seas), the optimized response is shown in Figure 39. The results of the four optimizations are 


shown in Table 17. 
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Figure 39. Simulation with sliding mode integral control in sea state three (head seas) 
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Optimized Values 2 I 











KT 
0 0 0 0 
0.0468/0.0464 0.0461/0.0453 0.025/0.025 0.0436/0.0553 
H/F (10? pounds) 19.6/-0.1 19.6/0.0 20.0/0.0 AO 
Mean Depth (feet) 55.01 55.02 55.16 


RMS Error (feet) 0.2345 (50%) 0.713 (49%) 0.693 (98%) 1.47 (11%) 
Maximum Error (feet) 0.789 (69%) 1.90 (45%) 1.92 (95%) 4.15 (13%) 


Sliding surface -0.1090 + 0.43141 -0.2418 + 0.28711 -0.7834 + 0.32601 
eigenvalues -0.1090 - 0.43141 -0.2418 - 0.28711 0075540927601 


0 0 0 0 





-0.0500 -0.0481 -0.0446 


Table 17. Optimized sliding mode integral control law results and performance 
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D. CONCLUDING REMARKS 

A comparison between the quality of control achieved by sliding mode control could 
conclude that the sliding mode control was inferior to state feedback control. This comparison 
would, however, neglect the added benefits of sliding mode control. The robust character of 
sliding mode control with the ability to provide reliable control for submarine control given 
uncertain hydrodynamic coefficients has been demonstrated for the NPS autonomous 
underwater vehicle program (Hawkinson, 1990). 

The sliding mode optimizations did not substantially reduce the control chatter and 
attendant high actuation rates. Variations of the sliding mode boundary thickness did not 
alleviate the chatter. 

Table 18 gives a summary of the RMS error achieved by each of the sliding mode 
optimizations. For comparison, it also includes the full state feedback results from Chapter IV. 
Although these were larger than the corresponding full state feedback cases, the sliding mode 
control proved to be much more robust in response to step changes in commanded depth. The 
sliding surface eigenvalues exhibited much more damping than the corresponding cases of full 
state feedback control. Also, it seemed to provide a more realistic average pitch angle for 


periscope depth operations. 








Sea State / Direction 3/Head 
Control Scheme S. с ~ 


RUE soe 
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Full State integral 0.372 0.536 
Sliding mode with integral control 0.2345 ЛЕ 0.693 


Table 18. Optimized RMS error (feet) of sliding mode control and full state feedback 
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УЕ GRAPHICAL DISPLAY 


A. INTRODUCTION 

In conducting ship control at periscope depth, a submarine diving officer relies on a 
variety of indications, meters, and some verbal reports to maintain the ship within the required 
depth band. In addition to the displayed parameters, the ship's control party also has their 
inertial reference, or “the seat of the pants”. It is perhaps the inertial reference which 
differentiates between great ship’s control parties, and the merely adequate. 

A submarine diving officer must track status of many ship's systems in addition to ship 


control. The items which the diving officer must monitor include: 


eMast positions 

eProximity of any portion of the ship to broaching (sail, rudder, mast fairing) 
eWater depth (general terms) 

eShip’s relationship to the submerged operating envelope 

eTrim 

eSpeed 

eWater density 

eShip’s evolutions (trash disposal, ventilation, etc.) 


e Towed array, floating wire antenna 


In many cases, the tracking tool most used is the diving officer's mental picture. 
Unfortunately, the ability to keep a clear status on many issues varies with fatigue and among 
individuals. This chapter gives the current conditions of the interface between the diving 


officer and ship's control, and proposes a different display medium to improve operations. 


B. CURRENT DIVING OFFICER INTERFACE 

To maintain a complete status, the diving officer has few tools at his disposal. He must 
rely on looking around at several different panels to get mast status, soundings, and water 
density while supervising the planesman. If an unplanned event, for example broaching, 


occurs the only record for reconstruction is the memory of the operators. 
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The gauges and meters used for ship control, while designed with generally appropriate 
time constants, are not adaptable for a given circumstance. For the most part, the same 
indications are used for high speed transit, periscope depth, and tactical operations. Figure 40 
and l'igure 41 show some of the indications used on board the USS Nautilus GSN 571). 


Although Nautilus 1s now a museum, the design of submarine ship’s control panels has not 
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Figure 40. USS Nautilus planes position indications 
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Figure 41. USS Nautilus pitch angle indication 


The current system of ship status display is very reliable, with redundant indications for 
important items and some purely mechanical indicators. It falls short in the area of presenting an 
integrated status. It is the writers opinion that the display system degrades the performance of 


the ship’s control party. With some parameters not displayed and others not conditioned for the 
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ship’s operating mode, operation near the surface in heavy sca states is extremely difficult. 
Skilled operators rely upon the existing indications, as well as “the seat of their pants” to 
maintain depth. Even so, it is a solid accomplishment to keep from broaching during extended 
periscope depth operations. 

{уеп more complex are operations in shallow water. ‘lhe proximity to grounding 
complicates all aspects of ship control. ‘lhe diving officer must be constantly aware of the water 
beneath the keel available for casualty recovery. Because nonzero pitch angles will cause one end 


of the ship to be deeper than indicated, this must also be accounted for. 


C. PROPOSED DISPLAY 

To incorporate the desired indications in a single display, a radically different approach is 
taken. Rather than rely on meters and gauges for the state of the ship, a screen is used. Figure 
42 shows the proposed display. A crude version of this display was developed using the 
SIMULINK® Animation Toolbox® (Figure 43). 
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Figure 42. Proposed graphical display of submarine control status 
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Figure 43. SIMULINK® animation of depth, pitch angle, and planes angles 


By integrating the ship status into one display, numerous improvements can be 
realized. The relationship of the submarine to the bottom and the surface is clearly shown. 
With the bottom contour information from a database, the diving officer has a continuous 
sense of the ship's proximity to grounding. In addition as sounding data 1s obtained, it can be 
displayed. 

The use of a digital display paradigm allows the display to be modified to support 
different operating modes. Because of the relationship between ship control and safety, the 
settings would be chosen based on a commanding officer approved doctrine. This would 
allow the operators to adjust the display system to best fit needs, and adapt it to new 
circumstances or missions. Alerts and alarms could readily added as the situation warranted. 

To assist the diving officer in maintaining status on the wave forces, several bar graphs 
were added to show the net force that the ship’s angle and planes were applying at a given 


time. These quantities would be filtered to provide a relevant average. Provided the averaging 
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interval was appropriate, this would queue the diving officer to order trim changes in response 
to changing environmental conditions. 

The periscope video in the bottom left hand corner would provide critical feedback for 
the dive. This would make the scope’s position relative to the surface apparent (another 
indication of depth), and allow the diving officer to be somewhat aware of the tactical 
situation. A close or new contact would prompt the diving officer to review mast exposure, 
which is also on the same display. 

Safety of shallow water operations would be enhanced by presenting a clear picture of 
the ship’s relationship to the bottom. During evasive action, the ship’s control party and the 
Officer of the Deck would be working with common knowledge of available water beneath the 
keel, and the contour ahead of the ship. 

Ship’s status could be recorded, to allow playback for the reconstruction of unplanned 


events. Figure 44 shows a possible data architecture to support the display. 
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Figure 44. Graphical display data paths 


D. CONCLUDING REMARKS 

The integration of the pertinent ship parameters in one display should yield dramatic 
improvements in submarine periscope depth operations. The diving officers improved 
awareness should reduce fatigue levels, allow for slightly lower speeds for a given sea state 
(reducing mast feather), and enable a much more complete environmental picture for the ship's 
control party. This awareness should increase the confidence of the ship's control party during 
demanding shallow water operations, reduce the likelihood of grounding or broaching, and 


provide an improved level of support for the Officer of the Deck. 
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УП. CONCLUSIONS AND RECOMMENDATIONS 


А. CONCLUSIONS 

As additional states were added to the control laws, the level achieved by the control 
law optimizations generally improved. For the full state feedback control cases, the depth 
control exceeded, in the authors experience, what is achievable by manual control. 

Ship’s control parties use more information about the state of the ship than is available 
from the explicit indications. Also, the success of the disturbance feedforward control suggests 
that averaged net force and moment would be of value to the ship’s control party. 

Sliding mode control provided well damped dynamics, with a robust behavior in 
response to changes in commanded depth. Although sliding mode control did not achieve the 
very low RMS errors of full state feedback, the gains which it employed were smaller and more 


realistic. 


B. RECOMMENDATIONS 

Although very good depthkeeping was achieved with state feedback control, the 
optimization schemes used resulted in control surface chatter, and very high control surface 
rates. To improve the quality of the model and provide more realistic planesman action 


several features could be added to the control laws and optimization routines: 


e Investigate other sea states, speeds and sea directions 

e Incorporate control surface rate limits 

e Include control surface chatter in the optimization objective functions 
e Use of Kalman filtering to provide state estimation and filtering 


e Investigate the use of depth rate for feedback control in place of heave 


The application of a new display system to an operating submarine is a major 
undertaking. Recommended steps to find a new manual submarine depth control paradigm 


arc: 


• Application of system identification techniques to submarine operating data to 


Investigate the nature of the human control 


ІІІ)! 


e Trials of a display onboard an appropriate submarine and or a submarine dive trainer 
e Use of recorded submarine Operating data to provide for “instant replay” training of 
ship's control personnel 
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APPENDIX A 


All computer code and SIMULINK® models used in this thesis are available from 
Professor Fotis Papoulias, Naval Postgraduate School. The computer programs and 


SIMULINK models used were: 


Programs: 

SUBOFF.M!- iniualizes SUBOFF hydrodynamic coefficients 

AXNL.M!- performs nonlinear Ax--E[qq;wq] calculation 

SMLOR!.M- determines MIMO sliding mode control law using Utkin's method (LOR) 
WF INI.M- reads wave data files, processes 

WFORCE.M- calculates wave forces for a given depth and time 

SB. INI.M- initializes model variables for MIMO vertical plane submarine model 

SBI INI.M- initializes model variables for MIMO vertical plane submarine model with integral 
depth control 

SB_SM.M- calculates MIMO SM control law from Q matrix 

SB_SS.M- calculates MIMO SM control law from sliding surface 

SB_SMFF.M- determines the feed forward matrix for a given sliding mode control law 
SB_PD.M-state feedback control law 

SB_PDFF.M- determines the feed forward matrix for a state feedback control law 
OBJ2.M - Objective function for pitch / depth feedback opumizations 

OPT2A.M- Optumization program for OBJ2.M and sea state three (head seas) 
OPT2B.M- Optumization program for OBJ2.M and sea state three (beam seas) 
OPT2C.M- Optimization program for OBJ2.M and sea state four (head seas) 
OPT2D.M- Optimization program for OBJ2.M and sea state four (beam seas) 
OBJ2ff.M - Objective function for pitch / depth feedback with disturbance feedforward 
optimizations 

OPT2FFA.M- Optimization program for OBJ2FF.M and sea state three (head seas) 
OPT2FFB.M- Optimization program for OBJ2FF.M and sea state three (beam seas) 
OPT2FFC.M- Opumization program for OBJ2FF.M and sea state four (head seas) 


! Given after list of programs 
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OPT2ZFFD.M- Optimization program for OBJ2FF.M and sea state four (beam seas) 
OBJ2I.M - Objective function for pitch / depth feedback with integral depth control 
optimizations 

OPT2IA.M- Opümization program for OBJ2I.M and sea state three (head seas) 
OPT2IB.M- Optimization program for OBJ21.M and sea state three (beam seas) 
OPT2IC.M- Optimization program for OBJ21.M and sea state four (head seas) 
OPT2ID.M- Optimization program for OBJ21.M and sea state four (beam seas) 
OBJ3.M - Objective function for full state partial distribution feedback optimizations 
OPT3A.M- Optimization program for OBJ3.M and sea state three (head seas) 
OPT3B.M- Optimization program for OBJ3.M and sea state three (beam seas) 
OPT3C.M- Optimization program for OBJ3.M and sea state four (head seas) 
OPT3D.M- Optimization program for OBJ3.M and sea state four (beam seas) 
OBJ3FF.M - Objective function for full state feedback with disturbance feedforward 
optimizations 

OPT3FFA.M- Optimization program for OBJ3FF.M and sea state three (head seas) 
OPT3FFB.M- Optimization program for OBJ3FF.M and sea state three (beam seas) 
OPT3FFC.M- Optimization program for OBJ3FF.M and sea state four (head seas) 
OPT3FFD.M- Optumization program for OBJ3FF.M and sea state four (beam seas) 
OBJ3I.M - Objective function for full state feedback with integral depth control optimizations 
OPT3IA.M- Optimization program for OBJSI.M and sea state three (head seas) 
OPT3IB.M- Optimization program for OBJ3I.M and sea state three (beam seas) 
OPT3IC.M- Optimization program for OBJ31.M and sea state four (head seas) 
OPT3ID.M- Optimization program for OBJ3I.M and sea state four (beam seas) 
OBJ3.M - Objective function for full state feedback optimizations 

OPTA4A.M- Optimization program for OBJ4.M and sea state three (head seas) 
OPT4B.M- Optimization program for OBJ4.M and sea state three (beam seas) 
OPT4C.M- Optimization program for OBJ4.M and sea state four (head seas) 
OPT4D.M- Optimization program for OBJ4.M and sea state four (beam seas) 
OBJ4FF.M - Objective function for full state partial distribution feedback with disturbance 
feedforward optimizations 

OPT4FFA.M- Optimization program for OBJ4FF.M and sea state three (head seas) 
OPT4FFB.M- Optimization program for OBJ4FF.M and sea state three (beam seas) 
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OPT4FFC.M- Opumization program for OBJ4FF.M and sea state four (head seas) 
OPT4FFD.M- Optimization program for OBJ4FF.M and sea state four (beam seas) 
OBJ4I.M - Objective function for full state partial distribution feedback with integral depth 
control optimizations 

OPTAIA.M- Optimization program for OBJAI.M and sea state three (head seas) 
OPTAIB.M- Optmization program for OBJ41.M and sea state three (beam seas) 
OPT4IC.M- Optimization program for OBJ4I.M and sea state four (head seas) 
OPT4ID.M- Optimization program for OBJ41.M and sea state four (beam seas) 
OBJ7.M - Objective function for sliding mode control optimizations 

OPT7A.M- Optimization program for OBJ7.M and sea state three (head seas) 
OPT7B.M- Optimization program for OBJ7.M and sea state three (beam seas) 
OPT7C.M- Optimization program for OBJ7.M and sea state four (head seas) 
OPT7D.M- Optimization program for OBJ7.M and sea state four (beam seas) 
OBJ7FF.M - Objective function for sliding mode control with disturbance feedforward 
optimizations 

OPT7FFA.M- Optimization program for OBJ7FF.M and sea state three (head seas) 
OPT7FFB.M- Optimization program for OBJ7FF.M and sea state three (beam seas) 
OPT7FFC.M- Optimization program for OBJ7FF.M and sea state four (head seas) 
OPT7FFD.M- Optimization program for OBJ7FF.M and sea state four (beam seas) 
OBJ7I.M - Objective function for sliding mode control with integral depth control 
optimizations 

OPT7IA.M- Optimization program for OBJ7I.M and sea state three (head seas) 
OPT7IB.M- Optimization program for OBJ7I.M and sea state three (beam seas) 
OPT7IC.M- Optimization program for OBJ7I.M and sea state four (head seas) 
OPT7ID.M- Optimization program for OBJ7I.M and sea state four (beam seas) 


OBJ7.M - Objective function for full state feedback optimizations 


Models: 
SMSW.M- MIMO sliding mode control submarine control model, with wave forces and trim 
SMSWFF.M- Same as SMSW.M, with disturbance feedforward 


SMT.M- used for determining steady state response, does not return the x state 


PDSW.M- MIMO state feedback control submarine control model, with wave forces and trim 


ТАРА 


PDSWFF.M- Same as PDSW.M, with disturbance feedforward 


PDT.M- used for determining steady state response, does not return the x state 
Wave Force Data files: 


For the first order motions / second order forces, each displacement / force 15 given in 
terms of a complex number. 


Wave spectral data files 
[wave frequency (rad/sec),S(@) (), amplitude (feet), phase (rad)] 


DESC A Rl, 
SE Grebe 
SPEC CHEI 
SPEC DIXI 


First order motions 
[0 (rad/sec), Wencounter (rad/sec), surge, sway, heavy, pitch, yaw, roll] 


LMOT_A.TXT 
OAMI 
EMOT C TXT 
LMOT_D.TXT 


Second order forces 


[01 (rad/sec), (02 (rad/sec), Fs, Fy, Fz, Mx, Му, М, | 


PORSSRTAT 
FORC#B FXT 
POR -CIXI 
ХЕ 


118 


SUBOFF.M 


function [A,B,MMI,Emat, D,WS]=suboff(U); 


% SDV hydrodynamic data, TM 231-78 page 6 


ff=2; Yofudge factor to make planes authority realisuc for 300 ft sub 


g—32.2; Vo 
L=300; % feet 


p=1.94; Yo slug/ft^3; 


m 2:0:015296*0:5*p5*T55; 
I —0.001084*0.5*p*L/^5; 
La = 0.00108*0.5*р*Г^5 ; 
mxg  =-0.127467E-4; 


W=m*g; 
zgb=1; % feet 
xgb=0; ИО ЕЕ 


Mqdot  —-0.000860*0.5*p*L/^5; 
Mwdot  —-0.000561*0.5*p*L^4; 
Mq . —-0.003702*0.5*p*L^4*U; 
Mw | —0.010324*0.5*p*L^5*U; 


Mds ——-ff*0.002409*0.5*p*L^3*U^2; 
Mdb =-Mds/4; 


Zqdot —-0.000633*0.5*p*L^4; 
Zwdot -—-0.014529*0.5*p*L^3; 


Та  =-0.007545*0.5*p*L*3*U; 
Zw  =-0.013910*0.5*p*L"2*U; 


Zds  -—-ff*0.005603*0.5*p*L^2*U^2; 
а —Zds/2: 


% define mass matrix, Compute mass matrix Inverse 


ШУ 


mm=eye(4); 
mm(1,1)=m-Zwdot; 
mm(1,2)—-Zqdot; 
mm(2,1)=-Mwdot; 
mm(2,2)=ly-Mgdot; 


A=zeros(4,4); 
A(1,D Zw; 
ШЕП Ыс; 
A(2,1)=Mw; 
A(2,2)=Mg; 
A(2,3)=-zgb*W; 
А(3,2)=1; 
А(4,1)=1; 
А(4,3)=-0; 


В-| “аһ Zds;¡Mdb Mds;0 0;0 0]; 
MMI=inv(mm); 
A=MMI*A; 


B=MMI*B; 
Emat=MMI(1:2,1:2)*[m*zgb,0;0,-m*zgb]; 


> 


— 


% diameter and surface area calculation 
Јар аб Se 


WS-67.651*(L/13.9792)^2; 
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AXNL.M 

function ax=axnl(x) 

global Amat Emat 
n=max(size(x)); 

a= Amat; 
a(4,1)=Amat(4,1)*cos(x(3)); 
x(3)=sin(x(3)); 


ax=a*x; 


Bep —ax(t1:2)-- Emat*[x(2)^2;x(2)*x(1)]; 


P 


SMLOR.M 
function [k,s,ks,e] =smlqr(a,b,q) 


% [k,s,ks,e] =smlqr(a,b,q) determines the sliding 


% mode control law for the system 
% xdot=a*x+b*u where 

% u=-kx-kn*satsgn(sigma) and 

% sigma=s*x. q 1s a positive definite 
% symmetric weighting matrix 

% used in assigning error weights 


% (LQR) to the states. Uses Utkin's method 
% as detailed in Hawkinson pp10-17 


[n,m]=size(b); 


% do transformation if required 
1£ norm(b(m+1:n,1:m))>eps*0.5 


[t,b1]=qr(b); 
b1=b1(1:m,1:m); 
с=с; 

else 
t=eye(n); 
b1=b(1:m,1:m); 


end 


q-t*q*t; 
а=(*а*г; 


411-а(1:п,1:т); 
q12=q(1:m,m+1:n); 
q21=q(m+1:n,1:m); 
q22=q(m+1:n,m+1:n); 
all=a(1:m,1:m); 
al2=a(1:m,m+1:n); 
a21=a(m+1:n,1:m); 
a22=a(m+1:n,m+1:n); 


as—a22-a21*inv(q11)*q12; 
qs=q22-q21*inv(q11)*q12; 


kt=are(as,a21*inv(q11)*(2215 qs); 


c2=(inv(q11)*(q12+221*k0); 
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k=-finv(b1)*(al 1 +c2'*a21),inv(b1)*(al2+c2'*a22)]*t; 
s=rref([eye(m); с2]'*0)'; 
ks=-inv(s'*b); 


e=elg([a+[b1;zeros(n-m,m)]*k*r")); 


22 
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APPENDIX B 


MAPLE? Solutions 


Determination of MIMO state feedback control steady state 


Determination of MIMO sliding mode control steady state 
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This maple script determines the steady state response in the vertical plane 


using MIMO state feedback control. Constant disturbances are assumed and the 
linear equations of motion applied to find the steady state depth error. 1% 
is also assumed that both control deflections are less than maximum. 
EQ1 and EQ2 are the heave and pitch linear equations of motion 
во = ба Ро шта з“ебегаљ-р 11“ ага "ај тр]ј 2 “џи - 2 хга2 ) ка 

EOl :=all wu+al30+bl1] u^ dl 4 b12 u^ d2 4 Fd 
вое ца м ша енеке са а. 


EQ2 := 421 wu+a23 0+b21 u2 dl +b22 u* d2 + Md 


readlib(isolate): 


w :- theta*u; 
w:=0u 
The linear state feedback laws at steady state are: 


dl := K11*w+K13*theta+K14*z; 
Qh = Ri l Out KII О 2 
d2 := K21*w+K23*thetarKl4á*z; 


а2с- К210Ө8и-К230-К147 


After the application of the control laws, the equations of heavy and pitch 


become: 
EON 


all8u?4a1304 bli u^ (Kll Ou KI30- Kldz) * bl2u^ (K21 0u-- K2310 WOO 
+ Ка=0 
EO02=05 


a210u2+a230+b21 u2(K110u+K130+KI142)+b22 u2 (K2] 0 u + К2ЗӨ А 
+ Md = 0 


Remove z from ЕО1 апа ЕО2 
БКО ОЕ (сое ералаш) 20 210 


FI : bl2u* KI4 * bH u^ K14 
> F2:=coeff (collect (expand(EQ2),2),2,1); 


к 
EQ3:-simplify(EQ1-EQ2*F1/F2); | | 


EO3 :=la110u2 522 +а110и2 621 +611 и3 КИ 022 +611 и? К13 0 b22 + al3 0 b22 
ОЗОР 2230511 + Fd b22 + Fdb2] -Mdb12—Ma bl1 + blow Koen 
+b12u2 K230b21-=a210u2 b12-a210u2 b11—b21 u> K11 0 b12-b21 u2 КА 


— b22 u? K21 0b11—b22 u2 K23 0b11-—a230b12)/(b22 + b21) 
{solate (E03, theta) DS XT кс 


Ө = ( -Fd b22 — Fd b21] + Ма Ь12 + Ма Ь11 an u^ БАВИ u* р21 +0611 ТЕ K1igb22 
И К на р та ВО азр - b12u^ K21521] tb 2 u K EPS 
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УОЛУ ЛУ О УУЛУ A NAAA 


бө О NU) 


> thetass:-coeff(collect(EQ3,theta),theta,0)/coeff(collect(EQ3,theta),theta,1); 


Meise = (Fa boos И ВО Malba? = MEP Cal ue LIA ЖО 
ADA а AS EE O 
a KIS DT OT DID al РО КГ Ре КО» 


— b22 u) K21 b11 — b22 u? K23 b11 — a23 b12) 


Substitute thetass into EQl1 and solve for steady state z 
EUST -coeff(collect(EQl1,z),z,0)/coefE£(collect(EO1,2),z,1); 


ен азе риска КВ) р КОГО ЕКОЕ 


ШЕШ) = > > 
bl2u* KI4- bll u^ KI4 


ИИ СС :- 
> coeff(collect(temp,theta),theta,1)*thetass+coeff(collect(temp, theta), theta,0) 


os = Gin on и? (K11 u - KI3) - b12 u* (K21l u+ K23)+ all u2) 
(Fd b22 + Fd b21 - Md b12- Md b11) [Cor u2 K14+b11 uz K14) (all u2 622 


= ЛО eS DI ЕАН кто а ОЕ 
E A E рю 


ШО u? K13b12-b22 u> K21 b11-b22 u2 K23 b11-a23 b12)) 
Ed 


b12 u^ KI4 * bll u? KI4 


B the coefficients Cl and C2 such that zss=C1*Fd+C2*Ma 
Еріп Ғу(ссеҒҒ(со11есг(ехрапа(255),Ға),Ға,1)); 


ZEE Dn v2-922 4 20] — a2Pu? b12- a21 49 511 — a23 012 2813 221 
NE ii 2. G13 b22 РОБ и? К 622 +2 Би” КЗ 22 — 21? КАБ 
ВР К21 621 22 и” К23'21 — 622 и’ КТ Ь - 622 и" K23 DTI 
— ати? К13 12 +612 и? К23 522 +11 и3 К11Ь21+ Ь12 иЗ К21 622 + 511 и2 К13621) 


н? KIA (128 b10) lali uz 622 + all ue 621 4 ОШ ОО и? КЗ 622 
3 


y 


+ 


KOT b2] + b по В 
КІЗ ОЕТ R20 


22 + 013 b21=a23 DI +b12 и 
E О пои 


– а23 b12)) 
> C2:=simplify (coeff (collect (expand(zss) ,Md) ,Md,1)); š 


с2 = - (аза + Ki pipu KI tb иЗ К21 + bl2u? K23 att |(( 








no 





аи” 622 +а1 и? 621 +11 и? К11 Ь22 +11 и К13 b22 + а13 Ь22  a13 b21 — a23 bli 


И ООСОР ee | 


— b21 u? K13 b12 — b22 u? K21 b11 — 522 u? K23 b11 — a23 b12) u? K14) 


Check that zss=C1*Fd+C2*Md 
> eq6:=simplify(zss-Cl*Fd-C2*Md); 


едб:=0 
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This maple script determines the steady state response in the vertical plane 
with basic MIMO sliding mode control. Constant disturbances are assumed and 
the linear equations of motion applied to find the steady state depth error. 
fies also assumed that both control deflections аге lecs hanna: ino rm 
mooie neither sliding SsUrface 1s saturated, 
EQ1 and EQ2 are the heave and pitch linear equations of motion 
EEUU -(lall*w*ucal3*theta-«-bll*u^2*dl4blI2*u^2 *d295bEd^ 
EQI :zall wu*al130- bl u* dl & b12 u^ d2 4 Fd 


EE M — (a7 baw *u+a23 *theta+b21*u°2 *d1i+b22 *u“*2*d2) +Md> 


КО = а Ри usas 2d u* dl + b22 u* d2 + Ма 


res dalib(isolate): 


БЕТ :- theta*u; 
и := ди 


The linear state feedback laws at steady state are: 


> di := K11*w+K13*theta+Ksl1l*etal*sigmal/phil+Ks12*eta2*sigma2/phi2; 

Ks11 etal sigmal " Ks12 eta2 sigma2 
phil phi2 

K21*w+K23*theta+Ks21*etal*sigmal/phil+Ks22*eta2*sigma2/phi2; 

Ks21 etal sigmal " Ks22 eta2 sigma2 
phil phi2 


dl = hl) Ou + Kis 6+ 


es a2 : 


dup А ОРН 22 ВЕ 


По оста |. =51] "+ 513" ћеба+514% 2; 

sremag = У и MOS 
РО он - 52] “+ 523 “"Ећега+524% 5; 

урад? S2 DENIED а А 


After the application of the control laws, the equations of heavy and pitch 
become: 


Б КӨШ: 


а11 К^ ТУ Ори [ки ӨТЕ ое ое USOS E) 


phil 


И Ks12 еіа2 (521 о ы 2 (кого и+ кэз 


phi2 
Ks21 etal (S11 0 u- 8513 0-- SI4z) | Ks22 eta2 (821 O.u + 523 0 + 524 7) 
E ll — ml 9 
phil phi2 





== w — J yy2RR, VZ II⁄IWKWKWQWIZIIZIIÉCI.I:€IIAIIIIII@IIIIIITWUsYA I T. TJ TC'UqCT2—U U T 





РАНО; = 
Кз] 
a21 0 u? - a23 0 b21 [ga + T SH HO 
рт 


Ks12 е1а2 (521 а as а (котви+ козе 


phi2 
Ка ега ичь/ 3 915132) Ks22 ciaz А о) 
fe A | + Ма = 0 
phil phi2 


Е-Е на 2... A 
ИЕ | = сое Collect (expand EOT) г),2,1); 
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— bll u* Ksl2 е1а2 524 , blu? Ksll eal $14. b12 u^ Ks22 eta2 $24 
phi2 phil phi2 
КР ЛЕВ 
phil 
>be. =COCTE (collect (expama(EhOZ) z), z, 1); 


. b21u? Ksl2 eta2 524 . b21 u? Ksll etal S14 ,222 u^ Ks22 eta2 S24 


T 


Г I + 
phi2 phil phi2 
252 uz K321 eta $4 
phil 


> E03'"=simpl1ty(E0l-=E02*E1/E2)5 
ВОЗЕ === E u^ КӨЛ" арала а сау к T K11 0 phil b22 Ks22 eta2 524 
— b11 u* Ks11 eta1 $13 0 b22 Ks22 eta2 S24 — b11 u) Ks12 eta2 $21 0 b22 Ks21 etal 514 
— b11 u* Ks12 eta2 S23 0 b22 Ks21 etal S14 — b11 u) K11 0 phi2 b22 Ks21 etal S14 
- b11 u? K13 0 phil b22 Ks22 eta2 S24 — b11 u^ K13 0 phi2 b22 Ks21 etal S14 
—a11 0 u? phil b21 Ks12 eta2 S24 — all 0 u^ phi2 b21 Ks11 etal S14 


– а110 u^ phil b22 Ks22 eta2 5324 — a110 u^ phi2 b22 Ks21 etal 514 
—al30 phil b21 Ks12 eta2 $24 — a13 0 phi2 b21 Ksl11 etal $14 
— al13 0 phil Ь22 К522 eta2 824 — a13 0 phi2 b22 Ks21 etal S14 


- b12 и? К21 9 рїї1 Б21 Кз12 еїа2 524 — b12 u? K21 0 phi2 b21 Ks11 etal S14 
DUE u^ K23 ВВ ПОД Матео 2 = 012 n“ К23 Өрйі2 b21 Ksll eta Ы 
— b12 u? Ks21 etal S11 Ө 221 К512 еіа2 524 – Б12 u^ Ks21 etal S13 0 b21 Ks12 eta2 S24 


= Ри и Ко22061а2 52 КӨ РАСК etaT Sd —b17 u^ Ks22 eta2 S23 0 b21 Ks11 etal S14 
– Fd phil b21 Ks12 eta2 $24 — Fd phi2 b21 Ks11 etal S14 — Fd phil b22 Ks22 eta2 824 


— Fd phi2 b22 Ks21 etal 514 + а21 0 u? phil b11 Ks12 eta2 $24 
t a21 0 u? phi2 b11 Ksll etal S14 a21 0 u^ phil b12 Ks22 eta2 S24 


t a21 0 u? phi2 b12 Ks21 etal S14 -- a23 0 phil b11 Ks12 eta2 S24 
+ a23 8 phi2 b11 Ks11 etal S14 + a23 O phil b12 Ks22 eta2 S24 


+ a23 0 phi2 b12 Ks21 etal S14 +621 u> K11 8 phil b12 Ks22 eta2 524 

552143 КІП Өрһі2 Ь12 К821 etal $14 & b21 u^ K130 phil b12 Ks22 eta2 524 

-b21u* K130 phi2 b12 Ks21 etal S14 + 621 иЗ Кз11 etal $11 0 b12 Ks22 eta2 $24 

+ 521 y? Ksll etal S130 b12 Ks22 eta2 S24 + 621 T Ks12 eta2 S21 0 b12 Ks21 etal S14 
& b21 u? Ks12 eta2 S230 b12 Ks21 etal SI4 4 b22 u? K21 0 phil b11 К512 eta2 S24 

t b22 u) K210 phi2 b11 Ks11 etal $14 b22 u^ K23 0 phil b11 Ks12 eta2 524 
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+ b22 u^ K230 РО РТК ета ТЕСЕ? и? КИНЕ ОУ ГРО БЕК КА 524 
+ 522 u? Ks21 etal S13 O0 bl11 Ks12 eta2 S24 —- D22 u^ Кора тета“ 32. ОРИ etal S14 
+ 522 и? Ks22 еіа2 523 9011 К511 ега! S14 + Md phil bl11 Ks12 eta2 $24 


+ Ма phi2 bll Ks11 etal $14 + Md phil b12 Ks22 eta2 S24 - Md phi2 b12 Ks21 etal $14)/( 
b21 Ks12 eta2 S24 phil  b21 Ks1l etal S14 phi2 + b22 Ks22 eta2 $24 phil 
+ 522 Ks21 etal S14 phi2) 
и omate(EO3,theta); 
9 = (Fd phil b21 KsI2 eta2 824 * Fd phi2 b21 Ksl11 etal S14 - Fd phil b22 Ks22 eta2 $24 
+ Fd phi2 b22 Ks21 etal S14 — Md phil b11 Ks12 eta2 S24 — Md phi2 b11 Ksll etal S14 


A О СОКО ОС у! 

-b11 u? Ks1] etal $11 b22 Ks22 eta2 $24 — b11 u^ K11 phil b22 Ks22 eta2 $24 

—bll s Ks12 eta2 521022 Ks2l etal 514- pII и? кола eta2 522 D22 KSO elal S14 
— all u2 phil b21 Ks12 eta2 $24 — b11 u? Ks11 etal S13 b22 Ks22 eta2 S24 

– БИ u3 K11 phi2 b22 Ks21 etal S14 — b11 u^ K13 phil b22 Ks22 eta2 S24 

— b11 u* K13 phi2 b22 Ks21 etal $14 — a13 phi2 b21 Ks11 etal S14 

—al1 u* phi2 b21 Ks11 etal S14 — all u^ phil b22 Ks22 eta2 S24 

—al1 u? phi2 b22 Ks21 etal S14 — a13 phil b21 Ks12 ега2 524 

- 512 и7 К23 рһ12 b21 Ks11 etal $14 — a13 phil b22 Ks22 eta2 S24 

— a13 phi2 b22 Ks21 etal S14 — b12 иЗ” K21 phil b21 Ks12 eta2 524 

— b12 u? K21 phi2 b21 Ks11 etal S14 — b12 u^ K23 phil b21 Ks12 eta2 S24 

- a21 u? phi2 b12 Ks21 etal S14 - a23 phil b11 Ks12 eta2 $24 

+ a23 phi2 b11 К511 еа! 514 — bI2 u? Ks21 etal $11 b21 Ksl2 eta2 524 

— 12 u2 Ks2] etal S13 b21 KsI2 eta2 824 — b12 u? Ks22 eta2 $21 b21 Ks11 etal S14 
— b12 u? Ks22 eta2 S23 b21 Ks11 etal 514 + а21 u^ phil b11 Ks12 eta2 S24 

- a21 u? phi2 b11 Ks11 etal $14 - a21 u^ phil b12 Ks22 eta2 524 

4 521 u? Ks12 eta2 S21 b12 Ks21 etal S14 -- b21 u^ Ks12 eta2 S23 b12 Ks21 etal S14 


+ 522 u? K21 phil b11 Ks12 eta2 524 + 522 75 K21 phi2 611 Ks11 etal S14 
+ a23 phil b12 Ks22 eta2 S824 + a23 phi2 b12 Ks21 etal S14 


Бори KI! phil b12 Ks22 eta2 S24 + 621 u’ KII phi2 b12 Ks21 etal S14 
4 b21 u? K13 phil b12 Ks22 eta2 $24 +b21 u* K13 phi2 b12 Ks21 eial S14 
В есть, Kosa a T Ро а 
+ b22 u? K23 phil b11 Ks12 eta2 S24 + b22 u? K23 phi2 KS etal S14 


ДЕЛІ 


3 
3 


ПОО ҮҮТ ҮҮ ОУ Т о. 


FOL? u eseectiazeszdebl 1 Ksi leetal S 1448022 u^ Ks22 eta2 $23 bl1 Ksil era $14) 


Find the steady state value of theta by setting EQ3=0 


> thetass:=simplify (coeff (collect (EQ3, theta) , theta, 0) /coeff (collect (EQ3, theta), 
нее lie) 


thetass := (Md phi2 b11 Ks11 etal S14+ Md phil b12 Ks22 eta2 S24 
+ Md phi2 b12 Ks21 etal S14 - Fd phil 621 Ks12 eta2 S24 + Md phil bI1 Ks12 eta2 S24 


_ Fd phil b22 Ks22 eta2 $24 — Fd phi2 b22 Ks21 etal S14 — Fd phi2 b2? Ks]! er ДУ, 
-b11 u? Ks1 etal S11 b22 Ks22 eta2 $24 — b11 u) K11 phil b22 Ks22 eta2 $24 

— b11 u? Ks12 eta2 S21 b22 Ks21 etal $14 — b11 u? Ks12 eta2 $23 b22 Ks21 etal S14 
RS phil b21 Ks12 eta2 S24 — bl1l u^ Kslletal S13 b22 Ks22 eta2 $24 

— b11 u? K11 phi2 b22 Ks21 etal $14 — b11 u? K13 phil b22 Ks22 eta2 $24 

— b11 u? K13 phi2 b22 Ks21 etal S14 —a13 phi2 b21 Ks11 etal S14 

— all u? phi2 b21 Ks11 etal S14 — all u^ phil b22 Ks22 eta2 S24 

— a11 u? phi2 b22 Ks21 etal S14 — a13 phil b21 Ks12 eta2 524 

— b12 u? K23 phi2 b21 Ks11 etal S14 — a13 phil b22 Ks22 eta2 $24 

— a13 phi2 b22 Ks21 etal $14 — b12 u? K21 phil b21 Ks12 eta2 S24 

— b12 u? K21 phi2 b21 Ks! etal SI4 — b12 u^ K23 phil b21 Ks12 eta2 S24 

- a21 u? phi2 b12 Ks21 etal $14 4 a23 phil b11 Ks12 eta2 S24 

+ 423 phi2 b11 Ks11 etal SI4 — bI2 u? Ks21 etal S11 b21 Ks12 eta2 S24 

— b12 u? Ks21 etal S13 b21 Ks12 eta2 S24 — b12 u? Ks22 eta2 S21 b21 Ks11 etal S14 
— b12 u? Ks22 eta2 S23 b21 Ksll etal S14 a21 u^ phil b11 Ks12 eta2 S24 

t a21 u? phi2 b11 Kslletal S14 a21 u? phil b12 Ks22 eta2 524 

& b21 u3 Ks12 eta2 S21 b12 Ks21 etal S14 4 b21 u* Ks12 eta2 $23 b12 Ks21 etal S14 


& b22 u? K21 phil b11 Ks12 eta2 S24 + b22 и? К21 phi2 b11 Ks11 etal S14 
+ a23 phil b12 Ks22 eta2 524  a23 phi2 b12 Ks21 etal $14 


+ Бо] из КГ phil b12 Ks22 eta2 S24  b21 u? K11 phi2 b12 Ks21 etal S14 

+ b21 u2 K13 phil b12 Ks22 eta2 S24 + Ь21 и? КІЗ phi2 b12 Ks21 etal S14 

+ b21u> KsI 1 etal S11 b12 Ks22 eta2 $24 - b21 u^ Ks11 etal $13 b12 Ks22 eta2 $24 
+ b22 u2 K23 phil b11 Ks12 eta2 524 + b22 u? K23 phi2 b11 Ks11 etal S14 

& b22 u) Ks21 etal S11 b11 Ks12 eta2 824 + b22 u? Ks21 etal $13 b11 Ks12 eta2 $24 
+ 122 u3 Ks22 eta2 S21 b11 Ks11 etal S14 ^ b22 u^ Ks22 eta2 $23 b11 Т N 514) 











Substitute thetass into EQ1 and solve for steady state z 


— — 





[02 


temp: =coett£tcallecr(EGL zz soe telde ы 


temp := (anı ди a13 0 


къта [kirous kso te SIL Out S130). КУ и (S27 0+ 250) 


plul phi2 
+ b12 u? È K210 u+ K3 op FELLS OUt SIZO) , K322 eta? (S218 u+ 5230) 
phil phi2 


AMA 


t D 2 
2 | Kslle al СУР ШШШ СЫР eta2 $24 | P | Ks21 eta] SI4  Ks22 ега2 524 
phil phi2 phil phi2 
S295 i- 
> coeff(collect (temp, theta) ,theta,1) *thetass+coeff(collect (temp, theta), theta, 0) 
7 Ë NU [ки Пи Ks11 etal (511 и+513) a Ks12 eta2 (521 u m 
phil phi2 
Ks2 са УЕ КБ22 ега2 (521 и =: 5280 D 
> AA P allu ( 
phil phi2 
Md phi2 b11 Ks1l etal S14 + Md phil b12 Ks22 eta2 524 + Md phi2 b12 Ks21 etal S14 
— Fd phil b2] Ks12 eta2 S24 + Md phil b11 Ks12 eta2 S524 — Fd phil b22 Ks22 eta2 524 


t b12 u* (кә! и+ К23 + 


— Fd phi2 b22 Ks21 etal S14 — Fd phi2 b21 Ks11 etal ғн) || 


Ш ‚2 | Ks11 etal 514 А Ks12 eta2 A. "TENE Кѕ21 еіа1 514 К522 еіа2 524 || 
phil phi2 

-b11 u? KsI1 etal S11 b22 Ks22 eta2 $24 — b11 u^ K11 phil b22 Ks22 eta2 $24 

—bll 9 Ks12 eta2 S21 b22 Ks21 etal SI4 — bll u* Ks12 е!а2 521022 4521 ва ton 

—a11 u? phil b21 Ks12 eta2 S24 — b11 u* Ks11 etal S13 b22 Ks22 eta2 524 

- b11 u? K11 phi2 b22 Ks21 etal S14 — b11 u^ K13 phil b22 Ks22 е1а2 524 

—b11 u* K13 phi2 b22 Ks21 etal S14 — a13 phi2 b21 Ks11 etal S14 

- all u2 phi2 b21 KsI] etal S14 — al1 u^ phil b22 Ks22 eta2 524 

—al1 u? phi2 b22 Ks21 etal S14 — a13 phil b21 Ks12 eta2 S24 

— b12 u? K23 phi2 b21 Ks11 etal S14 — al3 phil b22 Ks22 eta2 S24 

— a13 phi2 b22 Ks21 etal S14 — b12 u? K21 phil b21 Ks12 eta2 S24 

— b12 u? K21 phi2 b21 Ks11 etal SI4 — bI2 u* K23 phil b21 Ks12 eta2 524 

+ a21 u2 phi2 b12 Ks21 etal S14 + a23 phil b11 Ks12 eta2 $24 

+ a23 phi2 b11 Ks11 etal S14 — b12 u? Ks21 etal S11 b21 Ks12 eta2 524 

— b12 u? Ks21 etal S13 b21 Ks12 eta2 524-512 и? Кө22 eta2 S21 b21 Ks11 etal S14 


phil phi2 
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Б ООЗ Тата all ul pil bi? Ksi2 eta2 S24 
+ a21 u* phi2 b11 KsI] etal SI4 + а21 и? ри] b12 Ks22 eta2 S24 
Л УЛ О ОУУ УУ О. 


+ 522 u? K21 phil bI1 Ks12 eta2 524 + Ь22 иЗ K21 phi2 b11 Ks11 etal S14 
Tae Pil Dee KSZZ ela? 524-ра25 рін РЕЗЕ] стани“ 


t b21 и? K11 phil b12 Ks22 eta2 $24 + b21 u? K11 phi2 b12 Ks21 etal 514 
t b21 u* K13 phil b12 Ks22 eta2 S24 5521 u? K13 phi2 b12 Ks21 etal S14 
t b21 u? Ksl1l etal S11 bI2 Ks22 eta2 24 - b21 u^ Ksil eral S13 b12 Ks22 etad S24 
+ b22 u2 K23 phil b11 Ks12 eta2 S24 4 b22 и? K23 phi2 b11 Ks11 etal S14 


+ 522 ae КУЕ УМРЛА ПР? u? Ks21 etal SI3 b1l1 Ks12 eta ae 


2 


+ b22 We Ks22 eta2 S21 bll Ksll etal S14 + b22 u7 Ks22 eta2 S23 bll KST Eo n, 54) 


ға 
UECUIUHNIIITUOCEIUIENTUO | 
phil phi2 phil phi2 
Determine the све иселе" сл апа С2 SUCH that 2=6=6 | Ваа Иа 
> “Cl. = Б ву совебЕецео! Лесе емрапа 255). Еа). ва, ь)): 
C1: - (2 phi22 phil al3 522 Ks21 eral S4 2 phi2 pP a18 b21 KS12 ОКК, 
+ 2 phi2 phil? a13 b22 Ks22 eta2 524 4 2 phi2^ phil a13 b21 Ks11 etal 514 
+ phi2 phil? b11 u) K11 b21 Ks12 eta2 S24 +2 phi2 phil? b11 u) K11 b22 Ks22 eta2 S24 
* phi2 phil b11 u^ Ks11 etal S13 b21 Ks12 eta2 S24 
+ 2 phi2 phil b11 u^ Ks11 etal S13 b22 Ks22 eta2 524 
+ 2 phi2 phil b11 u? Ks11 etal S11 b22 Ks22 eta2 $24 
+2 phi2? phil b11 u? K13 b22 Ks21 etal S14 * phi2? phil b11 u^ K13 b21 Ks11 ега1 514 
+2 phi2 phil? b11 u? K13 b22 Ks22 eta2 524 + phi2^ phil b11 u? K11 b21 Ks11 etal 814 
+ phi2 phil b11 u? Ks11 etal 811 521 К512 ега2 524 
+ phi2 phil? b11 u? K13 b21 Ks12 eta2 S24 +2 phi2 
E A AI a T 
+ phi2 phil b11 u? Ks12 eta2 $21 b21 Ks11 etal 514 
+ phi2 phil b11 u? Ks12 eta2 $23 b21 Ks11 etal S14 
+ 2 phi2 phil b11 u^ Ks12 eta2 $23 b22 Ks21 etal S14 


+ 


> phi] b11 u? K11 b22 Ks21 etal S14 


+2 phi2 phil? b12 u? K21 b21 Ks12 eta2 $24 
+2 phi2 phil b12 u? Ks21 etal S11 b21 Ks12 eta2 S24 
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& phi2 phil? b12 u^ K23 b22 Ks22 eta2 524 + 2 phi2 phil^ b12 u^ K23 b21 Ks12 ега2 524 


4 2 phi2? phil b12 u? K21 b21 Ksl1 etal 514 * phi2^ phil b12 u? K21 b22 Ks21 etal S14 
+ phi2 phil? b12u? K21 b22 Ks22 eta2 S24-& 2 phi2? phi] bI2 u^ K23 b2d KsH etal S14 
+ phi2^ phil b12 u^ K23 b22 Ks21 etal S14 

+ phi2 phil b12 u* Ks21 etal S13 b22 Ks22 eta2 524 

+2 phi2 phil b12 u* Ks21 etal S13 b21 Ks12 eta2 524 

+ phi2 phil b12 u? Ks21 etal S11 b22 Ks22 eta2 524 

* phi2 phil b12 u? Ks22 eta2 S21 b22 Ks21 etal S14 

+2 phi2 phil b12 u? Ks22 eta2 S21 b21 Ks11 etal S14 

+ phi2 phil b12 u^ Ks22 eta2 523 b22 Ks21 etal S14 

2 phi2 phil b12 u^ Ks22 eta2 S23 b21 Ks11 etal S14 

$2 phi phil all u° b21 Ksl2 eta2 S24 +2 phi2 phil? all u® b22 Ks22 eta2 524 
T2mhi27 phil a11 u^ b22 Ks21 etal S14 +2 phi27 phil a11 u^ b21 Ksll etal S14 

X phi2? b11 u? Ks11 etal^ S11 b22 Ks21 S14 * phi2? b11 u? Ks11? etal? S11 b21 S14 
phil? b11 u? Ks12? eta2? $23 b21 824 4 phil? b11 u^ Ks12 eta2? 823 b22 Ks22 $24 
* phi2? b11 u? Ks11 eia1^ S13 b22 Ks21 S14 ^ phi2? b11 u? Ks11^ etal^ $13 b21 S14 
+ phil® b11 u? Ks12^ eta2^ S21 b21 $24 * phil? b11 u? Ks12 eta2^ S21 b22 Ks22 $24 
+ phi2? b12 u2 Ks212 etal” $13 b22 S14 ^ phi2? b12 u^ Ks21 etal? S13 b21 Ks11 S14 
& phil? b12 u? Ks22 eta2? 821 b21 Ks12 824 4 phil^ b12 u? Ks227? eta2? $21 b22 824 
4 phi2? b12 u? Ks21? eta? $11 b22 S14 ^ phi2? b12 u? Ks21 etal? $11 b21 Ks11 S14 
t phil? b12 u? Ks22 eta2^ S23 b21 Ks12 824 * phil^ b12 u^ Ks22^ eta2^ $23 b22 524 
— phi2? phil a21 u? b12 Ks21 etal S14 — phi2 phil? a23 b11 Ks12 eta2 $24 

— phi2? phil a23 b11 Ks11 etal S14 — phi2 phil^ a21 u^ b11 KsI2 eta2 S24 

— phi2? phil a21 u^ b11 Ks11 etal S14 — phi2 phil? a21 u^ b12 Ks22 eta2 $24 

— phi2 phil b21 u? Ks12 eta2 S21 b12 Ks21 etal S14 

— phi2 phil b21 u* Ks12 eta2 $23 b12 Ks21 etal S14 

— phi2 phil? b22 u? K21 b11 Ks12 eta2 824 — phi2? phil b22 u? K21 b11 Ks11 etal S14 
— phi2 phil? a23 b12 Ks22 е1а2 524 – phi2^ phil a23 b12 Ks21 etal S14 

— phi2 phil? b21 u? K11 b12 Ks22 eta2 S24 — phi2* phil b21 u? K11 b12 Ks21 etal S14 
— phi2 phil? b21 u? K13 b12 Ks22 eta2 $24 — phi2^ phil b21 u^ K13 b12 Ks21 etal S14 


2 
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— phi2 phil b21 u? KsI] etal S11 b12 Ks22 eta2 524 

— phi2 phil b21 u^ Ks11 etal S13 b12 Ks22 eta2 S24 

— phi2 phil? b22 u^ K23 b11 Ks12 eta2 S24 — phi2? phil b22 u? K23 b11 Ks11 etal S14 
— phi2 phil b22 u? Ks21 ега1 S11 b11 Ks12 eta2 524 

— phi2 phil b22 u* Ks21 etal S13 b11 Ks12 eta2 S24 

— phi2 phil b22 u? Ks22 eta2 S21 b11 Ks11 etal S14 

=p DA n^ Ks22 eta2 S23 b11 Ks11 etal $14) (и? (phil b11 Ks12 eta2 S24 


+ phi2 b11 Ks} I etal S14 + phit 12 Ks22! 009804 CpRIIDT? o P СД, 

bll u?’ КӨП etal 911 622 Ks22 6102 $24 РТТ nno am T 

— b11 u? Ks12 eta2 821 b22 Ks21 etal S14 — b11 u^ Ks12 eta2 $23 b22 Ks21 RESB 
-all u” phil b21 Ks12 eta2 $24 — b11 u^ Ks11 etal $13 b22 Ks22 е1а2 524 

- b11 u? K11 phi2 622 Ks21 etal S14 — b11 u* K13 phil b22 Ks22 eta2 $24 

- b11 u* K13 phi2 b22 Ks21 etal S14— al3 phi2 b21 Ks11 etal 514 

ЛА y? phi2 b21 KslIl etal S14 —all y? phil b22 Ks22 eta2 524 

— all u” phi2 b22 Ks21 etal 514 — a13 phil b21 Ks12 е1а2 524 

- b12 u2 K23 phi2 b21 Ks11 etal $14 —a13 phil b22 Ks22 eta2 S24 

—al3 phi2 b22 Ks21 etal S14 — b12 u” K21 phil b21 Ks12 eta2 S24 

- b12 u? K21 phi2 b21 Ks11 ета1 514 — b12 u? K23 phil b21 Ks12 eta2 S24 

+ a21 u* phi2 b12 Ks21 etal S14 + a23 phil b11 Ks12 eta2 $24 

+ a23 phi2 b11 Ks11 etal S14 — b12 u? Ks21 etal S11 b21 Ks12 eta2 S24 

— b12 u? Ks21 etal S13 b21 Ks12 eta2 S24 — b12 u? Ks22 eta2 S21 b21 Ks11 etal S14 
— b12 u? Ks22 eta2 S23 b21 Ks11 etal S14 * a21 u^ phil b11 Ks12 eta2 S24 

+ a21 u2 phi2 b11 Ks11 etal S14 * a21 u^ phil b12 Ks22 eta2 524 

+ b21 u? Ks12 eta2 $21 b12 Ks21 etal SI4 + b21 u2 KsI2 eta2 S23'612'Ks21 SN 


1 o из K21 phil b11 Ks12 е!а2 524 + 622 и? K21 phi2 b11 Ks11 etal S14 
+ a23 phil b12 Ks22 eta2 $24 + a23 phi2 b12 Ks21 etal S14 


t b21 u? K11 phil b12 Ks22 eta2 524 * b21 u” K1] phi2 b12 Ks21 etal S14 
+ b21 u2 K13 phil b12 Ks22 eta2 524 * b21 u^ K13 phi2 b12 Ks21 eta1 $14 
& b21 u9 Ks1l etal S11 b12 Ks22 eta2 S24 4 b21 u* Ks11 etal S13 b12 Ks22 eta2 $24 
& b22 u? K23 phil b11 Ks12 eta2 S24 & b22 u? K23 phi2 b11 Ks11 etal S14 
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+ 522 из Ks2] etal Sil ofl Ksi2 etaz P 5 o и? КАД ааа ОТКОР Od 


+ b22 u” Ks22 eta2 $21 b11 Ks11 etal 814  b22 u^ Ks22 eta2 $23 b11 Ksl1 etal 514)) 
> €C2:- simplify(coeff(collect(expand(zss),Md),Md,1)); с 


о = (кај! etal phi2 b11 и? $11 + Ksil etal phi2 bl] u? S13 + etal phi2 u^ $13 b12 Ks21 
+ etal phi2 u? S11 b12 Ks21  phi2 b11 phil u^ K11 * phi2 b11 phil u? K13 + phi2 phil a13 
+ phi2 phil b12 u^ K23 * phi2 phil a11 u?  phi2 phil b12 w^ K21 

Y b1 phil u? Ks12 eta2 $21 & b11 phil u? Ks12 eta2 $23 + phil b12 u2 Ks22 eta2 523 


+ phil b12 ТЕ Ks22 eta2 s21)/((o11 и? Кеа КОКА Еа 74 


- b11 u? K11 phil b22 Ks22 eta2 S24 — b11 u^ Ks12 eta2 S21 b22 Ks21 etal S14 
— b11 u* Ks12 eta2 S23 b22 Ks21 etal S14 — all u? phil b21 Ks12 eta2 $24 

— b11 u* Ks! eral S13 b22 Ks22 eta2 $24 — b11 u^ K11 phi2 b22 Ks21 etal S14 
— b11 u* K13 phil b22 Ks22 eta2 S24 — b11 u^ K13 phi2 b22 Ks21 etal S14 
Кро D2] Ksil eral Sl4—all и“ phi2 b21 Ks11 etal S14 

– а11 u* phil b22 Ks22 eta2 S24 — a11 u* phi2 b22 Ks21 etal S14 


— a13 phil b21 Ks12 era2 S24 — b12 u* K23 phi2 b21 Ks11 etal S14 
— al13 phil b22 Ks22 eta2 $24 — a13 phi2 b22 K$21 etal 514 


— b12 u? K21 phil b21 Ks12 eta2 $24 — b12 u? K21 phi2 Ь21 К511 etal S14 


— 512 u^ K?23 phil b21 Ks12 eta2 824 + а21 а phi2 b12 Ks21 etal S14 
+ a23 phil b11] Ks12 eta2 S24 * a23 phi2 b11 Ks11 etal 514 


— b12 u? Ks21 eal S11 b21 Ks12 eta2 $24 — b12 u^ Ks21 etal S13 b21 Ks12 eta2 S24 
— b12 po Ks22 eta2 S21 b21 Ksll etal S14 — b12 u* КУ? еа 57 а Ба У 
+ а21 и“ рћи b11 Ks12 eta2 $24 + а21 и? phi2 b11 Ks11 etal S14 

1 a21 u? phil b12 Ks22 eta2 S24  b21 u? Ks12 eta2 S21 b12 Ks21 etal S14 

22 Ks]2'eta2 $23 b12 Ks21 eta! SI4-& b22 u? K21 phil b11 Ks12 eta2 S24 

+ b22 u> K21 phi2 b11 Ks11 etal S14 + a23 phil b12 Ks22 eta2 S24 

+ a23 phi2 b12 Ks21 etal S14  b21 w? K11 phil b12 Ks22 eta2 S24 

i202 u^ КІ! рһ12 12 К521 ега! 514 + Б21 иг КІЗ phil b12 Ks22 eta2 524 

+ b21 u2 K13 phi2 b12 Ks21 etal S14 * b21 u? Ks11 etal S11 b12 Ks22 eta2 $24 

4 521 u? Ks1 etal S13 b12 Ks22 eta2 S24 - b22 u^ K23 phil b11 Ks12 ега2 524 

+ b22 u2 K23 phi2 b11 Ks11 etal $14 22 и Ks21 etal S11 b11 Ks12 eta2 S24 
ИО ООО enl ST DIT KsI2ea2 824 39022 w^ K322 eta SD Rea 
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Check that zss-Cl*Fd-«C2*Md 
> едо = 5 ]ifyvy zss Cl+pa- Co n 


ego Ü) 
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